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ABSTRACT 

Context. The Galactic Cold Cores project has carried out Herschel photometric observations of 116 fields where the Planck survey 
has found signs of cold dust emission. The fields contain sources in different environments and different phases of star formation. 
Previous studies have revealed variations in their dust submillimetre opacity. 

Aims. The aim is to measure the value of dust opacity spectral index and to understand its variations spatially and with respect to 
other parameters, such as temperature, column density, and Galactic location. 

Methods. The dust opacity spectral index (5 and the dust colour temperature T are derived using Herschel and Planck data. The 
relation between p and T is examined for the whole sample and inside individual fields. 

Results. Based on IRAS and Planck data, the fields are characterised by a median colour temperature of 16.1 K and a median opacity 
spectral index of P- 1.84. The values are not correlated with Galactic longitude. We observe a clear T-p anti-correlation. In Herschel 
observations, constrained at lower resolution by Planck data, the variations follow the column density structure and ySpiR can rise to 
~ 2.2 in individual clumps. The highest values are found in starless clumps. The Planck 111 GHz band shows a systematic excess that 
is not restricted to cold clumps and is thus consistent with a general flattening of the dust emission spectrum at millimetre wavelengths. 
When fitted separately below and above 700//m, the median spectral index values are Pyir ~ 1.91 and P(mm) ~ 1.66. 

Conclusions. The spectral index changes as a function of column density and wavelength. The comparison of different data sets and 
the examination of possible error sources show that our results are robust. However, p variations are partly masked by temperature 
gradients and the changes in the intrinsic grain properties may be even greater. 

Key words. ISM: clouds - Infrared: ISM - Submillimetre: ISM - dust, extinction - Stars: formation - Stars: protostars 


1. Introduction 

The all-sky survey of the Planck satellite ( [Tauber et al.||2010| 
enabled the identification of cold interstellar dust clouds on a 
Galactic scale. The high sensitivity and an angular resolution 
of better than 5' allowed the detection and classification of a 
large number of cold and compact sources that, based on the 
low temperature alone, are likely to be associated with prestellar 

* Planck {httpillwww.esa.intIPlanck) is a project of the European 
Space Agency - ESA - with instruments provided by two scientific con¬ 
sortia funded by ESA member states (in particular the lead countries: 
France and Italy) with contributions from NASA (USA), and telescope 
reflectors provided in a collaboration between ESA and a scientific con¬ 
sortium led and funded by Denmark. 

** Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and with im¬ 
portant participation from NASA. 

Table 3 is available in electronic form at the CDS via anony¬ 
mous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u- 
strasbg.fr/cgi-bin/qcat?J/A-hA/ 


phases of the star formation process. Analysis of Planck data led 
to the creat ion of the Cold Clump Catalogue o f Planck Objects 
(C3PO, see |Planck Collaboration XXIII|201l] ), with more than 
10000 sources. The low colour temperatures (mostly T ^ 14 K) 
indicate that the fields have high column density structures that 
are only partially resolved by the ~ 5' Planck beam. Many of the 
clumps are likely to contain gravitationally bound cores. 


The Herschel Open Time Key Programme Galactic Cold 
Cores carried out dust continuum emission observations of 116 
fields that were selected based on the Planck C3PO catalogue. 
The fields were mapped with Herschel PACS and SPIRE instru¬ 


ments ( jPilbratt et al.||2010] Poglitsch et al.||201^ [Griffin et~ar 
2010[ ) at wavelengths of 100-500 yum. Herschel makes it possi 
ble to study the Planck clumps and their internal structure in de 


tail ([Juvela et al.|20T^[Montillaud et al.|2015, in press[[Rlvera- 


Ingraham et al.[[2015, submitted[ ). The mapped fields are typ¬ 
ically ~ 40' in size, enabling investigation of dust properties 
on larger scales and, in particular, comparison of dust proper¬ 
ties between the dense and cold regions and their lower den- 
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sity environment. First results have been presented in Planck 


[Collaboration XXIl 

[I[(201 l[);[Planck Collaboration XXIIj([201 1[), ( 

and in Juvela et al. 

POTO MTT[2012|i (Papers I, II, and III, re- i 

spectively). [Montil 

aud et al. (2015, in press]) present an analysis 

of submillimetre c 

umps and star formation in these Herschel i 


fi elds. Other stud i es hav e been carried out in high latitude fields 
Malinen et al.| ( [2014| ), |Rivera-Ingraham et al.| ( |2015, submit-| 


in 


ted|), and Ristorcelli et al. (in preparation). 


In this paper we concentrate on the dust opacity spectral in¬ 
dex, p. Our previous studies indicated significant variations in 
the dust submillimetre opacity, k, and the average opacity was 
found t o be more than tw ice the typical value found in diffuse 
clouds ( [Juvela et'ST 20I5| ). In the densest clumps - if nearby and 
hence well resolved - k showed additional increase by more than 
a factor of two. We therefore might expect to see some related 
changes in the spectral index for these sources. 

The observed spectral index or emission spectral index, p, 
is derived from the spectral energy distribution (SED) of dust 
emission assuming that dust is optically thin, and its emission 
can be modelled with a modified blackbody. In such a model p 
describes the asymptotic behaviour of the dust opacity at long 
wavelengths. According to physical models describing the opti¬ 
cal response of matter with light, such as the Lorentz model, at 
far-infrared (FIR) and submillimetre wavelengths yS is indepen¬ 
dent of temperature and may take values i n the range 1- 2 de- 
pending on dust composition and structure ( jBohren & Huffman] 
|I998| ). However, FIR/submillimetre observations have shown 
that (i) p values depend on the astrophysical environment, (ii) 
that p and the dust temperature are anti-correlated, and (Hi) that 
P varies with wavelength such that the dust SED flattens for 
A > 500-800 yum. Th e yg- T anti-correlation w as first observed 
by jPupac et al.| ( 2003] ) and jPesert et al.| ( |2008) ) from PRONAOS 
and ARCHEOPS observations of various regions of the interstel¬ 
lar medium (ISM) in the FIR and submillimetre-to-millimetre 
ranges respectively. Many studies based on Herschel and Planck 
data find that p and Tdust are anti-correlated (jRodon et al.|2010 


Ven eziani et al.|20T0l [Paradis et al.|2()T0l|Etxaluze et al. 2011 

Planck Collaboration et al.|2011 1. These studies have derived p 

values as low as ~ 1.5 in warm regions, while p values greater 
than 2 have been obtained for cold environments. 

A flattening of the dust SED at longer wavelengths, first 
observed by [Reach et al.[ ( [1995[ ), is also widely detected by 
Herschel and Planck in various Galactic and extragalactic envi¬ 
ronments ([Gordo n et al. [20101 [Galliano et al. [20 lT| [Paradis et al. 
2012[ [Planck Collaboration XVII[[2011t [Planck (Collaboration 
XIX|2011[ ). In the literature, this is sometimes referred to as the 
sub-mm or FIR excess. Rather than a single value of p across the 
entire FIR-to-millimetre regime, one may thus define two spec¬ 
tral indices, ySpiR at v > 353 GHz and p^m at v < 353 GHz. This 
approach was adopted by Planck studies of ISM dust at high 
galactic latitude ([Planck Collaboration Int. XVn|2Q14 ) and in 
the Galactic plan ^Planck Collaboration Int. XIV|2014 ), which 
found that ySpiR is greater than p^m by ~ 0.15 and ~ 0.3 dex, 
respectively. 

The observed p - T anti-correlation could stem from the in¬ 
trinsic properties of dust grain material, more specifically from 
its amorphous structure. Laboratory studies of interstellar dust 
analogues show that the dust opacity increases with tempera¬ 
ture, that the spectral index increases when the temperature de¬ 
creases, and that a single spectral index is generally not suf¬ 
ficient to describe the experimental data ([Agladze et al][1996[ 


Mennella et al.|1998|[Boudet et al.|2005 ; [Coupeaud et al.|201 1[ ). 
This temperature-dependent behaviour is related to low energy 
transitions that are associated with the amorphous structure of 


the material. It can be modelled with the TLS model ([Meny 
erar][2007l ), which has also been shown to successfully repro¬ 


duce observations of ISM dust emission at FIR/submillimetre 
wavelengths ( Paradis et al. [2011 2014 ). Variations in the dust 
spectral index could also reflect the dust evolution in the various 
astrophysical environments, e.g. grain growth and coagulation in 
dense cloud s ( [Steinacker et al. 12010} [Pagani et al.|2010[ [Kohler 
et al. [2012|) or evolution of the carbonaceous dust component 


( [Jones et al.|2013[ ). It is probable that several processes occur in 
the ISM and the interpretation of the T-p anti-correlation is there¬ 
fore not straightforward. Finally, it is also important to note that 
P and Tdust are partially degenerate. Estimates of p are affected 
by noise in the data, temperature variations along the line-of- 
sight and by the method used to fit the observations ([Shetty et al. 


2009b|al [Juvela & Ysard[[2012a 

Malinen et al.[[2011l [Juvela & 

Ysard|2()12bl Ysard et al.|20121 

Juvela et al.|20131 [Pagani ^1- 


2015). The robustness and interpretation of any observed anti¬ 
correlation between p and T^ust requires the accurate determina¬ 
tion and proper treatment of these biases. 

The structure of the paper is the following. The observations 
are described in Sect.j^ The main results are presented in Sect.[^ 
including estimates of the average value of p and its dependence 
on dust temperature. The results are discussed in Sect. We 
summarise our conclusions in Sect. O 


2. Observations and methods 

2.1. Herschel data 

The target fields for Herschel observations were selected using 
the information that was available from Planck observations and 
ancillary data at that time. The goal was to cover a represen¬ 
tative set of Planck clumps with respect to Galactic longitude 
and latitude, estimated dust colour temperature, and clump mass. 
The selection of the Herschel fields is described in[ 


et al. 


( 2012| ) and an overview of all the maps is given in 
( [2015, in press] ) and [Juvela et al.[ ( [2Q15| ). 


uvela et all 


Montillaud 


he sample 

o 

does not include sources in the Galactic plane, \b\ < 1 , which 


is covered by the Hi-GAL programme ( [Molinari et al. 2010). 
Similarly, regions included in other Herschel key programmes 
like the Gould Be lt survey ( [Andre et al.[[2010| ) and HOBYS 
( [Motte et al.|2010| ) were avoided. 

Our observations cover 116 fields. The SPIRE maps corre¬ 
spond to wavelengths 250yum, 350yum and 500jam. They have 
an average size of ~ 1800arcmin^. PACS maps were obtained 
at 100/im and 160yum. They are smaller, with an average size 
of ~ 660arcmin^. The fields are liste d in Table [H Tht Herschel 
observation numbers can be found in [Montillaud et al.[ ( [2015, in] 
[press] ). 

The SPIRE observations were reduced with the Herschel 
Interactive Processing Environment HIPE v.12.0, using the of¬ 
ficial pipeline with the iterative destriper and the extended emis¬ 
sion calibration options. The maps were produced with the naive 
map-making routine. The PACS data at 100yum and 160yum were 
processed with HIPE v. 12.0 up to Level 1 and the final maps 
were produced with Scanamorphos v.23 ( [Rousse l|2013p . In or¬ 
der of increasing wavelength, the resolution is approximately 7", 
12", 18", 25", and 37" for the five bands. The raw and pipeline 
reduced data are available via the Herschel Science Archive. 
The accuracy of the absolute calibration of the SPIRE obser¬ 
vations is expected to be better than 7% and the relative calibra- 
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tion between the SPIRE bands better than 2^^ For PACS, when 
combined with SPIRE data, we adopt an error estimate of 10%. 
The value is consistent with the differences between PACS and 
Spitzer MIPS measurements of extended emission]^ 

Herschel observations were convolved to the resolution of 
the 500yum data using the latest information about the beam 
shape^ The effective point spread function (PSF) depends on 
the source spectrum. For spectra characteristic of our fields (tem¬ 
peratures -10-20 K, spectral index -1.5-2.5) the variation of the 
PSF shape at a given distance from its centre is at a level of 
- 1 %. The net effect on the final surface brightness maps is 
smaller, because of radial averaging and because the signal is 
never produced by a single point source. The difference with the 
Gaussian approximation used in the previous papers (e.g. |Juvela] 
et al.|201^|Montillaud et al.|2015, in press] ) is more significant 


Planck (IRIS and the 857 GHz Planck data were convolved to 


this resolution; see below Sec t.]2.21). The detailed procedure o f 
the comparison is described in Montillaud et al.l ( 2015, in pre^ . 
The Planck and IRIS measurements are interpolated to Herschel 
wavelengths using modified black body curves with (3 = 1.8. 
Finear least squares fits between Herschel and these reference 
data provide estimates of the Herschel surface br ightness zero 
points and their statistical uncertainty. Unlike in jJuvela et~ST 



Fig. 1. Correlations between Herschel and reference data in 
the field G300.86-9.00. PACS 160yum is plotted against AKARI 
140yum data (left frame, V resolution) and SPIRE 250yum data 
against interpolated IRIS and Planck data (right frame, 5.0' res¬ 
olution). 


(see IGriffin et al. 112013] ). We use beams calculated for a modi¬ 
fied blackbody spectrum with fixed parameters T = 15.0 K and 
P = 2.0. These correspond to the expected parameter values for 
cold and compact sources but the beam shape does not vary sig¬ 
nificantly for the range of T and p values found in our fields. 
We constructed circular symmetric convolution kernels that, to¬ 
gether with the original PSF, result in an effective PSF identical 
(to within a fraction of one per cent) to that of the 500 yum band. 
After convolving data to the common resolution of the 500 yum 
band, all maps were resampled onto the 14" pixel grid of the 
original 500yum maps. This provides adequate sampling of the 
- 37" beam. For part of the analysis in this paper, the data were 
further convolved to 1.0' resolution. 

To determine the absolute zero point of the intensity scale 
of Hersc hel S PIRE data, we c ompared the data with Planck 
(see S ect. 12.2] ) and IRIS 100yum (Miville-Deschenes & Fagache] 
2005] ) observations. This was done at the -4.85' resolution of 


with AKARI 140 yum all-sky maps ( ]Doi et al.]2015] ) with a resolu¬ 
tion of - 1' and with the surface brightness scale tied to DIRBE 
(resolution - 40'). The comparison includes colour corrections 
and extrapolation between 140yum and 160yum. Both are calcu¬ 
lated from the same modified blackbody spectra as above, using 
yS = 1.8 fits to SPIRE channels (see below). The 160yum gain 
calibration is not changed and we continue to assume a rela¬ 
tive uncertainty of 10%. The comparison of Herschel diXid AKARI 
intensities (averaged over the Herschel coverage) involves very 
small statistical errors. The uncertainty of the 160yum zero point 
is thus dominated by the uncertainty of the AKARI zero point, the 
relative calibration, and the wavelength extrapolation. We sub¬ 
tracted from AKARI data a cosmic infrared background (CIRB) 
level of 0.95MJysr“^ but the un certainty of this comp onent is 
estimated at only -0.2MJysr“^ ( Matsuura et al.]]201l| ). We as- 


( 12012] ), the results are used to correct not only the zero point 
but also the gain calibration. At 250 yum the estimated zero point 
errors are -IMJysr"^ or less. The formal uncertainties of the 
slopes (gain correction) are below 1% and typically -0.5%. We 
continue to assume an uncertainty of 2% for the calibration of 
the SPIRE channels. Figure shows an example of the quality 
of the correlations between PACS 160yum and AKARI 140yum 
and between SPIRE 250 yum and interpolated IRIS and Planck 
data. Note that at the lowest surface brightness the PACS error 
estimates are already dominated by the multiplicative calibra¬ 
tion uncertainty, which is 4MJy sr“^ for a surface brightness of 
40MJy sr“^ (see below). 

For PACS, the linear correlations are often less well defined. 
In our previous papers (e.g. ] Juvela et~ar]12011| 12012] ), the zero 
points were determined from direct comparison of the average 
surface brightness values of the Herschel and interpolated IRAS 
and Planck data. Because of the similarity of the wavelength, 
here we use zero points that are obtained from a comparison 


^ SPIRE Observer’s manual, 
httpiUherschel esac. esa. int/Documentation. shtml 
^ http://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalibrationWeb 
^ http://herschel.esac.esa.int/twiki/bin/view/Public/ 
SpirePhotometerBeamProfileAnalysis 


sume that the overall zero point accuracy of 160 yum is not better 
than -1 MJy sr“\ the statistical uncertainty of the zero points in 
the previously used correlations with Planck and IRIS. The rea¬ 
son for using Herschel 160yum data instead of relying directly 
on AKARI 140yum band is that, in addition to higher spatial res¬ 
olution, the map quality of Herschel data is better, especially at 
small scales. 

The zero point and gain corrections discussed above were 
calculated iteratively together with the determination of colour 
corrections. In the case of SPIRE, the corrections also take the 
dependence between source spectrum and beam solid angles into 
account. The colour corrections were calculated for modified 
blackbody spectra with yS = 1.8 (the expected average value over 
the fields), using colour temperatures that were calculated pixel- 
by-pixel using SPIRE data and the same fixed value of p. At 
T = 15.0 K the difference between p = 1.8 and p = 2.2 colour 
corrections would be -1.0% for all channels 160yum-500yum 
(the lOOyum band is not used in this paper). Part of that error 
would further cancel out when a change in p is compensated by 
an opposite change in the fitted temperature. Even in the 160yum 
band and at a low temperature of 12 K, a difference of Ap = 0.4 
in the colour correction translates to a change of less than 1 % in 
surface brightness. Thus, colour corrections with a fixed value 
of p do not bias the p values that we derive in our subsequent 
analysis. 

The SPIRE error estimates mainly consist of the assumed 
uncertainty of gain calibration that is tied to Planck via least 
squares fits that are performed separately for each field and band. 
Thus, these include a significant statistical component, in ad¬ 
dition to the uncertainty of absolute calibration of Planck (and 
IRIS). In the 250yum band, additional uncertainty results from 
interpolation between lOOyum and 350yum. The PACS uncertain- 
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ties consist of the uncertainty of the original Herschel gain cal¬ 
ibration, of zero point uncertainties (resulting from comparison 
with AKARI), and of mapping artefacts that can be more impor¬ 
tant than for SPIRE. The correlations between PACS and inter¬ 
polated IRAS and Planck data (not used in this paper) also sug¬ 
gest that PACS errors are not dominated by systematic errors. 
In SPIRE bands systematic errors should be a larger fractional 
component because of the smaller statistical errors. 

2.2. Planck and IRIS data 

Planck data from the 857 GHz, 545 GHz, 353 GHz, and 217 GHz 
bands are used to study the dust spectrum at lower resolution, 
to make use of the longer wavelength coverage of Planck. The 
857 GHz and 545 GHz channels were already used above to es¬ 
timate zero points for the Herschel surface brightness data. The 
IRIS 100/im data suffer from the problems mentioned earlier 
(see Sects. PI and |2.3| ), the unknown contribution from stochas¬ 
tically heated very small grains and the sensitivity to line-of- 
sight temperature variations. However, because the combination 
of Planck and IRIS has already been studied over the whole sky 
(e.g. [Planck Collaboration XII||2014| ), it is useful to investigate 
them also in the selected Herschel fields. 

Th e Planck data are described in I Planck Collaboration VIII I 
( |2Q14| ). We use the data from the official release that is avail¬ 
able via the Planck Legacy Arc/z/v^ and is described by the 
Planck Explanatory Supplemenj^ We use zodiacal light sub¬ 
tracted maps where the estimated cosmic infrared background 
has also already been subtracted. The cosmic microwave back¬ 
ground (CMB) anisotropies were subtracted using the CMB 
maps, which have been calculated with the SMICA algorithm 
( [Planck Collaboration XII|20T4| ) and are available in the Planck 
Legacy Archive. To set the Galacti c zero levels, the surface 
brightness values listed in Table 5 of [Planck Collaboration Vllf 
( [20I4[ ) were subtracted. These do not contain the cosmic in¬ 
frared background and, for consistency, an estimated CIRB level 
of 0.80MJy sr“^ was also subtracted from 100//m IRIS data 
( [Hauser et al.| 19981 [Matsuura et aI.|2Ql 1[ ). 

Several transitions of the CO molecule fall within the Planck 
bands, producing signal above the pure dust emission. The tran¬ 
sition / = 2-1 is inside the 217 GHz band and the transition 
J = 3-2 inside the 353 GHz band. The higher transitions are 
likely to be too weak to be a concern because of the steep rise of 
the dust spectrum towards shorter wavelengths. In cold clouds 
and regions of low mass star formation, t ypical peak intensities 
of J = 3-2 CO lines are ~10Kkms“^ (e.g. San Jose-Garcia et ^ 
2013| . This corresponds to -0.5 MJy sr“^ in the 353 GHz band, 
but dilution by the ~ 5' beam is likely to make the contribution 
smaller (CO peak values are measured at a resolution of - 1' or 
less). In the most prominent massive star-forming regions, the 
line intensities can be higher by a factor of ten and CO con¬ 
tributions could be locally significant in the 353 GHz band. At 
217 GHz, the line intensities are of similar order but the con¬ 
version factors from line intensity to Pla nck surface brightness 
is lower b y more than a factor of two ( [Planck Collaboration] 
XIII[[2014] ). However, because the dust signal decreases from 
353 GHz to 217 GHz by a factor of -5, CO contamination is 
more important in the 217 GHz band and, if not taken into ac¬ 
count, could bias the {3 estimates downwards. 

We do not have direct multitransition CO observations of our 
fields. One exception is LDN 183 (G6.03-r36.73), where the av¬ 


erage line ratio rA(2-l)/rA(l-0) over the central 3' x 6' area is 
-0.54. A second example is cloud LDN 1642 (G210.9 0-36.55) 
which has b een observed with the SEST telescope by [Russeil 


et al.[ ( [2003] ). SEST data gives an average ratio of line areas of 
W{2 - 1)/W(1 - 0) = 0.61. At a resolution of 4.5', one reaches a 
value of 0.72 at the location of the column density peak. 

Planck data themselves have been used to estimate, with 
three different methods, the intensity of the first CO transitions. 
The methods used to derive the C O estimates are described in 
Planck Collaboration XIII ( [2014[ ). Planck Type 1 CO maps are 


calculated based on the differences in the spectral response of 
bolometer pairs. This makes it possible to estimate the CO emis¬ 
sion of the three first rotational transitions separately, although 
with relatively low signal-to-noise ratios (SN). Type 2 estimates 
of CO(l-O) and CO(2-l) line intensities result from a multi¬ 
channel fitting of the CO emission together with cosmic mi¬ 
crowave background (CMB), dust emission, and free-free emis¬ 
sion. Because of the use of lower frequency channels. Type 2 
maps only have a resolution of 15'. Einally, Type 3 maps are 
calculated assuming fixed CO line ratios over the whole sky, 
CO(2-1)/CO(1-0)=0.595 an CO(3-2)/CO(l-0)=0.297. This re¬ 
sults in maps with the best SN but, of course, no separate in¬ 
formation of the line ratios in any given pixel. The reason why 
we do not make the correction for CO emission directly using 
Type 1 and Type 2 maps is the low signal-to-noise ratio (SN) 
of Type 1 maps (at full Planck resolution) and the low spatial 
resolution of Type 2 maps. The component separations also as¬ 
sume a fixed dust spectrum (T = 17K, yS = 1.6), which may 
lead to errors in regions that are colder than the average ISM. 
According to Planck Collaboration XIII ( 2014[ ), Type 1 maps 
give at large scales line ratios rA(2-l)/rA(l-0)=0.4 and 
2)/rA(l-0)=0.2. However, concentrating on CO bright regions, 
the same paper found average line ratios of 0.595 and 0.297. 

Eigure shows the line ratios we obtained in our fields for 
data averaged over the SPIRE coverage. Type 1 data give values 
0.51 and 0.25 for the (2-l)/(l-0) and (3-2)/(l-0) ratios. Using 
Type 2 data (with lower resolution but higher SN) the (2-l)/(l- 
0) ratio is 0.53. The values are thus in line with the values quoted 
above for large area averages and, for (2-l)/(l-0), close to the 
value observed in LDN 183. The values are somewhat lower than 


in CO brig ht regions in general (see above [Planck Collaboration 
XIIl|[2014[ ), which may be related to the fact that our fields are 


^ http://pla.esac.esa.int/pla/aio/planckProducts.html 
^ http://www.sciops.esa.int/wikiSI/planckpla 


selected based on their low (dust) temperature. In Pig.we have 
divided each field to two parts according to its median column 
density (estimated using SPIRE data and a constant value of yS = 
1.8). There are no clear differences in the line ratios between 
the two column density samples. However, Type 1 (3-2)/(l-0) 
ratio seems to increase in the high column density sample at the 
highest IT(l-O) values. A similar but weaker trend might exist in 
the (2-l)/(l-0) ratios of Type 1 data but is not seen in the Type 
2 data. 

In this paper we correct the 345 GHz and 217 GHz data for 
CO contamination using the CO(l-O) emission of the Planck 
Type 3 CO maps and as a default assuming line ratios Ta(2- 
l)/rA(l-0)=0.5 and Ta(2-2)/Ta(1-0)=0.3. The conversion fac¬ 
tor between line intensity and its contribution to surface bright¬ 
ness are given in Planck Collaboration XIII ( 2014[ ). The ratio 
Ta(3-2)/Ta(1-0)=0.3 is above the average observed ratio (see 
Fig. 13 but is representative of the ratios at the highest col¬ 
umn densities. The correction in the 353 GHz band is typically 
-0.1 MJy sr“\ at most a few per cent of the total signal. The se¬ 
lected ratio rA(2-l)/rA(l-0)=0.5 is consistent with the data in 
Pig.|^ The CO correction in the 217 GHz band can be more than 
10 %. 
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Fig. 2. CO line ratios in our fields based on Planck Type 1 and 
Type 2 CO maps. Type 1 ratios (2-l)/(l-0) and (3-2)/(l-0) are 
plotted with blue and red circles and Type 2 ratios (2-l)/(l-0) 
with green squares. The open and the filled symbols are averages 
over pixels below and above median column density, respec¬ 
tively. The black lines show least squares fits to data including 
both low and high column densities. The shaded areas indicate 
the range of line ratios considered in this paper. 


In |PIanck Collaboration XIII| ( |20I4|) the Planck Ty pe 3 
^^CO(l-O) maps were compared with Dame et al. (2001) data 
of Taurus, Orion, and Polaris. The ratios between the Planck 
and Dame et al. line intensities were found to be slightly above 
one, with an average value of 1.085. In case of Type 1 maps, 
the / = 2-1 estimates were compared with ground-based data 
and with estimates derived from COBE FIRAS, both compar¬ 
isons showing agreement within -10%. Another comparison of 
J = 2-1 observations is now possible u sing the large-sca l e map s 
of Orion A and B clouds presented in |Nishimura et al. ( |2015| ). 
Orion is not representative of our sample (selected based on the 
low dust temperature of the sources) but this has little effect on 
the comparison of the CO measurements. Over the area pre¬ 
sented in Fig. 1 of|Nishimura et al.|(|2015|), the average ratio be¬ 


tween the Planck Type 2 CO estimates and the 1.85 m Nobeyama 
radio telescope data is slightly below one, 0.91+0.11. The ratio 
does not show systematic variations as a function of column den¬ 
sity or dust temperature, although the value is lower, -0.6, at the 
location of Ori KF. The mean line value is compatible with the 
estimated uncertainty of -10 % of the intensity calibr ation of the 
Nobeyama data. Figure 11 in |Nishimura et al.| ( |2015| ) also shows 
that, in spite of large column density changes, the spatial varia¬ 
tions of the (2-l)/(l-0) line ratio are relatively small. In particu¬ 
lar, the ratio remains relatively constant -0.8 over most of Orion 
B. This suggests that the ratio is affected more by the general 
radiation field of the region than by the local column density. 

The comparison with ground-based data shows that, within 
an accuracy of -10%, the combination of Planck Type 3 and 
Type 2 data give adequate predictions of the average contamina- 


tionbyCO J = 1-0 and J = 2-1 lines (see [Planck Collaboration 
XIII|2014] ). Type 3 data directly trace the CO emission at 6' res 


olution. Thus, the main uncertainty in the local CO correction of 
the 217 GHz band results from possible variations of line ratios. 
The examples of FDN 183 and FDN 1642 show that even within 
very dense cold clumps the (2-l)/(l-0) ratio does not necessarily 
rise much above the mean value of -0.5. This in spite of the fact 
that for very high column densities the line ratio should tend to¬ 


wards one. Higher line ratios (^1) are possible in hot cores where 
the gas kinetic temperature is affected by local star-formation. In 
Orion, because of the high-mass stars, the (2-l)/(l-0) ratio is 
in many places close to one or even above one (see Fig. 11 in 
Nishimura et al.|20T5] ). Comparison with dust temperature maps 


(based on Herschel archive SPIRE maps and an assumption of 
P = 1.8) shows that the dust temperature of those regions is 
-20 K or above. At Herschel resolution, only 19 of our fields 
contain pixels with T > 20K. These amount to -3% of the to¬ 
tal mapped area. Temperatures T > 20 K are found mainly in 
regions of low column densities and some of the high values re¬ 
sult from higher uncertainties at low intensities and near map 
borders. Nevertheless, there can be some denser regions with 
both strong CO emission and high (2-l)/(l-0) line ratios where 
our default CO correction would be clearly too small. The inten¬ 
sity and temperature maps (see Appendices) can be consulted to 
identify regions with potential for such errors. 

If CO correction is underestimated, it can affect the derived 
dust spectral index and, in particular, estimates for p variations 
as a function of wavelength. We examine later the effect of ap- 
plyingCO corrections that are twice the default correction (see 
Sect. |4.4| ). This corresponds to an assumption of line ratios (2- 
1)/(1-0)=1.0 and (3-2)/(l-0)=0.6. These represent a firm upper 
limit for the average CO contamination. This also is an upper 
limit for cold clumps where the CO excitation temperature is be¬ 
low -15 K. Farger line ratios are still possible in regions of very 
intense radiation field. The shaded regions in Fig. [^illustrate the 
difference between the two assumptions of line ratios. 

The 217 GHz band could also contain a non-negligible 
amount of free-free emission. To estimate and to subtract this 
contribution, we used the low frequency foreground component 
model (available in the Planck Legacy Arc/t/v^ and its predic¬ 
tions at 217 GHz. The model consists of an estimate for the low 
frequency spectral index (Healpix maps with NSIDE=256 and 
angular resolution of 40') and higher resolution estimates of in¬ 
tensity (Healpix maps with NSIDE=2048 and full Planck reso¬ 
lution). These result from component separation where the other 
components are cosmic microwave background, CO emission, 
and thermal dust emission ( jPlanck Collaboration XII||2Q14| ). 
The component separation cannot take the local variations of 
the dust spectral index, dust temperature, relative intensity of 
synchrotron and free-free components, and the contribution of 
anomalous microwave emission fully into account. However, it 
is a useful indication of the level of non-dust contributions and 
its subtraction should reduce whatever small bias this emission 
might cause. For our fields, this component ranges from 0.02% 
to 2.3% of the total signal of the 217 GHz channel. The median 
value is 0.1% (for averages over the maps). 

For the instrumental noise, we use the error estimates in¬ 
cluded in Planck frequency maps. The accuracy of absolute cal¬ 
ibration is estimated to be 10% for the 857 GHz and 545 GHz 
channels and -2% or better at low frequencies. The error es¬ 
timates associated with the CO and low frequency corrections 
were set conservatively to 50% of the correction. For the noise 
caused by cos mic infrared background (GIB), we use the esti¬ 
mates given by [Planck Collaboration XI ( [2014[ ). All error com¬ 
ponents were added together in quadrature. At 857 GHz and 
545 GHz the error budget is dominated by calibration uncertain¬ 
ties while at 217 GHz the uncertainty of the CO correction is 
often the largest term. The instrumental noise and GIB fiuctua- 
tions have significant contribution only at low column densities, 
in the most diffuse parts of the fields. 


http://pla.esac.esa.int/pla/aio/planckProducts.html 
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The data between 100/im and 353 GHz were fitted pixel-by- 
pixel using modified black body spectra with a fixed spectra in¬ 
dex of p = 1.8 and the data were colour corrected for this spec¬ 
trum. For the range of colour temperatures derived from SPIRE 
observations and for spectral index changes of +0.2, a change in 
the Planck colour corrections could change the surface bright¬ 
ness values by ~1%. Thus, the subsequent analysis is only very 
weakly dependent on the value of {3 assumed in the colour cor¬ 
rection. 

The Planck beams are characterised in IPlanck Collaboration! 
|Vn| ([20T4]). The FWHM of the effective beam is 4.63', 4.84', 
4.86', and 5.01' for the 857 GHz, 545 GHz, 353 GHz, and 
217 GHz bands, respectively. The beams have some ellipticity 
that varies across the sky. The ellipticity is particularly large 
in the 857 GHz band, on average e ~ 1.39. Our calculations 
are based on the assumption of symmetric beams. In the zero 
point determination of Herschel data, Planck 857 GHz maps 
were convolved to 4.85' assuming Gaussian beams. This convo¬ 
lution does little to mask the effects of beam ellipticity, which 
thus causes additional noise in the determination of Herschel 
zero points and becomes part of the estimated statistical error. 
When Planck data are used in modified blackbody fits, all data 
are convolved to a resolution of 5.0'. We made some calculations 
at 8.0' resolution, where the effects of beam asymmetries should 
be small, but did not observe any significant discrepancy with 
the 5.0' calculations. 


The IRIS versions of th e IRAS 100//m maps ( [Miville- 
Deschenes & Lagache||2005] ) were convolved from the origi 


nal resolution of 4.3' to 4.85' (zero point calculations) or 5.0' 
(other analysis), assuming circular symmetric Gaussian beams. 
The IRIS data are corrected for the expected contribution of very 
small grains (VSGs) which raises the 100yum intensity above 
that of the pure big grain (BG) emissi on. The corrections were 
derived with the DustEM dust model ( Compiegne et al.||2011| ) 
as a function of the BG temperature (see Juvela et al.||2()12 ). It 
accounts for -10% at the highest temperatures and rises up to 
-20% for the coldest fields. The uncertainty of this correction 
is large but mostly within the 15% uncertainty adopted for the 
lOOyum band. The VSG correction mainly decreases the colour 
temperature estimates and the effect on the spectral index is 
small. 


2.3. Calculation of dust opacity spectral index 

The spectral index is determined by fitting observed surface 
brightness values /y with a modified blackbody law 

Iy = IoX (By(T)IByo(T)) X (v/vof. (1) 

The fit involves three free parameters: the spectral index p, the 
colour temperature T, the intensity /q at a reference frequency vq 
and By is the Planck function. The fit is possible when observa¬ 
tions consist of at least three wavelengths. Observations at long 
wavelengths (beyond 200 pm) are needed to constrain the spec¬ 
tral index, while shorter wavelengths (A ^ 200yum) are better 
for determining colour temperature. The equation includes the 
assumption of optically thin emission. 

The interpretation of the estimated I3(T) dependence requires 
special care, because errors can produce artificial correlations 
( [Shetty et al.|2009at [Juvela & Ysard|20I2a| ). Even small errors 
in surface brightness data may lead to significantly higher T and 
lower p or vice versa. One must also note that T and p only 
characterise the shape of the observed spectrum and will dif¬ 
fer from the corresponding intrinsic dust parameters. The colour 


temperature overestimates the mass-averaged physical dust tem¬ 
perature and the observed spectral index is smaller than the av- 
erage opacity spectral index of the grains ([Shetty et al.||2009b[ 

[Malinen et al.|201l| [Juvela & Ysard|20I2b| ). These differences 

are related to a mixture of grain temperatures along the line-of- 
sight or generally within the beam. These can be caused by tem¬ 
perature gradients inside the clouds (for example, cold clumps 
with warmer envelopes or embedded hot sources) or to the pres¬ 
ence of different grain populations that have different tempera¬ 
tures. Because wavelengths near and shortward of the peak of 
the emission spectrum are more sensitive to temperature, the re¬ 
covered parameters will depend on the set of wavelengths used. 
At -lOOyum and below, the contribution of stochastically heated 
VSGs increases and the spectrum cannot be approximated as a 
modified blackbody due to the broader temperature distribution 
of VSGs. The PACS lOOyum data are therefore not used in this 


paper. 

The least squa res fi ts of Eq . employ the error estimates 
given in Sections ^ and 2^ The uncertainties were calcu¬ 
lated with direct Monte Carlo simulations or with Markov chain 
Monte Carlo (MCMC) methods. With MCMC one can calculate 
the full posterior probability distributions of /q, T, and yS, only 
using the relative probabilities of the fits performed with differ- 
ent combinations of the p arameter values (e.g., |Veneziani et ak] 
|20IQt [Juvela et^|2Q13| ). The calculation results in a chain of 
parameter values and we use their median as the final parame¬ 
ter estimates. In selected fits, the zero point uncertainties were 
included as part of the error estimates (combined as additional 
statistical error) or their influence was examined separately as a 
systematic error affecting each field individually. 


3. Results 

We use different combinations of observational data to seek ev¬ 
idence for spectral index variations between the fields and be¬ 
tween regions characterised by different temperature and column 
density. We utilise some internally consistent data sets (e.g., a 
single instrument) and compare the qualitative and quantitative 
results obtained with different and even independent data sets 
(e.g., Herschel data vs. Planck and IRAS data). We start by us¬ 
ing Eq. [^to fit Planck data in combination with IRAS lOOyum 
band. Next we fit Herschel observations of the three SPIRE 
bands, 250yum, 350 yum, and 500yum. We continue with four 
bands, adding the PACS 160yum band to SPIRE data. We fin¬ 
ish by using a combination of Herschel and Planck data. In these 
fits, dust temperature is determined mainly by Herschel observa¬ 
tions, Planck providing further constraints on the spectral index, 
albeit at a lower spatial resolution. The different cases are listed 
in Table The table lists the adopted relative errors to which 
other sources of uncertainty (CO, GIB, and low frequency model 
subtraction) are added. 

3.1. Results from Planck and IRIS data 

We use IRIS lOOyum and Planck data in the 857 GHz, 545 GHz, 
353 GHz, and 217 GHz bands to estimate colour temperatures 
and spectral indices at scales of 5.0' and larger. Unlike in the 
case of Herschel data, we need to include the short wavelength 
IRIS band to constrain dust temperature. This is expected to be 
more sensitive to warm dust and may include contribution from 
stochastically heated small dust grains. Th ese effects have been 
discussed extensivel y in the literature (e.g. [Shetty et al.[[2009b[ 

[Malinen et al.|201l1 Juvela & Ysard|2012b[ ). Radiative transfer 

models of Appendix ^ suggest that the inclusion of the lOOyum 
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Table 1. Combinations of observational data used. 


Section Data Resolution Error estimates 


3.1 

IRIS lOOyL/m, Planck 857-217 GHz 

5.0' 

15%, 10%, 10%, 2%, 2% 

3.2 

Herschel 250, 350, 500/am 

37" or V 

2% + zero point uncertainties, which have mean values of 
0.39, 0.17, and 0.07 MJy sr"^ for the three bands 

3.3 

Herschel 160, 250, 350, 500jam 

37" or V 

10%, 2%, 2%, 2% + a zero point uncertainty, which 
for the 160 fim band is assumed to be 1 MJy sr"^ 

3.5 

Herschel 250-500 jjm and 

Planck 857-353 GHz 

- 37" 

Herschel 3.5%, 

Planck 10% (857-545 GHz) and 2% (353-217 GHz) 



Fig. 3. Colour temperature (upper frame) and spectral index values (lower frame) for the 116 fields. For clarity, alternating field 
names appear above and below the figure. The values are derived from IRIS and Planck data, averaged over a single FWHM=10' 
beam. Estimates are calculated with (red circles) and without (blue squares) the 217 GHz band. The black horizontal bars show 
values derived without the 217 GHz band, us ing surface brightness averaged o ver a larger beam with FWHM=30'. The fields are 
arranged in the order of increasing distance ( [Montillaud et al.||2015, in press] ). The five vertical lines indicate the first fields with 
estimated distances above 0, 200, 400, 800, and 2000 pc. The first eight fields do not have reliable distance estimates. 


band does not make the (3 determination much more sensitive to 
line-of-sight temperature mixing. For a cloud with ry = 10 mag, 
the combination of IRIS and the three Planck channels led to p 
values that were up to 0.07 units lower than for the combination 
of SPIRE and the same Planck bands, the total bias being -0.10 
units. 


We start by calculating one SED per field. In each field, 
we average surface brightness with a Gaussian beam with 
FWHM=10.0', located at the centre of the Herschel SPIRE cov¬ 
erage. Because the Herschel maps are typically 30-40' in diam¬ 
eter, the values are representative of most of the field and espe¬ 
cially of the higher column density regions usually found at its 
centre. SEDs were fitted with modified blackbody functions and 
the results are shown in Fig.[^ This figure includes estimates cal¬ 
culated with and without the 217 GHz band. The error bars have 
been estimated with Monte Carlo sim ulati ons, using the surface 
brightness uncertainties listed in Sect\22 


When the lowest frequency is extended from 353 GHz 
to 217 GHz, the mean spectral index is lower. The average 


change is Ap = -0.10 and the mean temperature is higher by 
AT = 0.52 K. The IRIS 100yum dat a hav e been corrected for 
the expected VSG emission (see Sect. |2.2| ). The inclusion of the 
VSG correction decreases temperatures by AT —0.5 K and in¬ 
creases spectral index values by Ap ~ 0.05. With the 217 GHz 
band included, the reduced values (normalised by the num¬ 
ber of degrees of freedom) are many times higher. Although 
the number depends on the relative error estimates of the bands, 
this already strongly suggests that a single modified black body 
curve does not provide a good fit over the full wavelength range. 

The horizontal black bars show values derived with IRIS 
and Planck bands down to 353 GHz band, using surface bright¬ 
ness measurements averaged over a beam with FWHM=30'. 
Compared to the 10' values (blue points), the temperatures are 
higher and spectral index values lower when estimated from the 
lower resolution data. This difference is expected because our 
fields are selected for the presence of cold dust and, therefore, 
typically have larger than average column densities. The emis¬ 
sion properties are thus clearly different already between the 
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scales of 10' and 30'. In |Planck Collaboration XXIlI| ( |2011| ), the 
temperatures observed towards individual Planck cold clumps 
were even lower (the distribution peaking below 14 K) and the 
spectral index values were higher, p ~ 2.1. However, that anal¬ 
ysis concerned individual clumps (sizes below ~10') and only 
the cold emission component that was obtained by subtracting 
an estimate of the extended warm background emission. 

Figure]^ shows the correlations between colour temperature 
and spectral index for the 10' beams, for values calculated with¬ 
out the 217 GHz band. The median values over all fields are 


T = 15.88 K and1.85. 

The distribution of (T, P) points in Fig. shows some neg¬ 
ative correlation between T and p that cannot be entirely ex¬ 
plained by measurement errors. The errors scatter {T, P) points 
along a direction that has a steeper negative slope than the over¬ 
all parameter distribution. In the surface brightness errors, with 
the possible exception of the 353 GHz and 217 GHz bands, the 
main component is the uncertainty of the absolute calibration. 
The band-to-band errors should be smaller, because calibration 
errors are likely to be correlated between the fields. To illustrate 
the potential importance of this point, we recalculated the error 
regions using 50% smaller uncertainties for the surface bright¬ 
ness measurements (thick contours in Fig.|^. 

In the figure, we overplot the relation p{T) = (0.4 -r 0.008 x 
T)~^ derived from PRONAOS observations of cold submillime¬ 
tre clumps ( Pupae et al.||2003[) and the corresp onding relation 


po 

P{T) = 11.5 X from [Desert et al.| ( |2008| ). Both relations 

cross the distribution of our data points. 

Figure shows an example of the fits performed pixel-by- 
pixel using the IRAS 100 jim and Planck 857 GHz, 545 GHz, and 
353 GHz data. The field is G300.86-9.00, part of the Musca fil¬ 
ament. In frame J, the spectral index is seen to increase towards 
the filament and especially towards the two dense clumps. The 
change is from p -1.75 in the background to a peak value of 
P = 2.03 in both clumps. In Fig. frame e shows the residuals 
between Planck 217 GHz data and the fit, and frame / shows the 
difference between the spectral indices estimated in frequency 
intervals 3000 GHz-353 GHz and 3000 GHz-217 GHz. These 
show that the spectrum is flatter at the longest wavelengths. In 
relative terms (frame /), the change of spectral index as a func¬ 
tion of wavelength appears to be larger in the diffuse medium. 
However, line-of-sight temperature variations tend to decrease 
the fitted value of p, thus reducing the 217 GHz excess. This may 
partly explain the spatial pattern in frame /. 

Figure shows correlations between T and p. These are 
calculated pixel-by-pixel, using IRIS and Planck maps (v > 
353 GHz) at 5' resolution. The figure shows a distribution similar 
to that of Fig.|^ This suggests that spatial averaging from 5' to 
10' scales has not significantly affected the appearance of Fig. 

In the second frame of Fig. we use only pixels with estimated 
column densities in the fourth quartile (25% of pixels with the 
highest column densities). To select the pixels, we use column 
densities that are calculated with the same four frequencies but 
employing a fixed value of the spectral index, p = 1.8. Figure 
shows that the distribution is similarly elongated at high column 
densities but has shifted towards lower temperatures and higher 
values of spectral index. The median values are T = 16.59 K and 
P = 1.80 for all pixels and T = 15.92 K and p = 1.86 for pixels 
with column density in the fourth quartile. 

Appendix shows the (T, P) fits for all 116 fields. A 
217 GHz excess (i.e. a flattening of the dust SEP at long wave¬ 
lengths) is present in most fields. In some fields (e.g. G202.23- 
3.38 and G253.71-rl.93 the situation is less clear and the residu¬ 
als are at least partly negative. There are an additional two fields 



Fig. 5. Fit of (r, P) in the field G300.86-9.00 using IRAS and 
Planck data. The frames a, c, and d show 250 yum surface bright¬ 
ness, colour temperature, and spectral index fitted to 3000 GHz 
(IRAS 100yum), 857 GHz, 545 GHz, and 353 GHz data, respec¬ 
tively. Frame b shows, for comparison, the higher resolution 
Herschel t( 250 yum) map. Frame e is the surface brightness resid¬ 
ual between the 217 GHz observation and the above fit. The last 
frame shows the difference between spectral indices derived with 
3000 GHz-353 GHz and 3000 GHz-217 GHz data, lip decreases 
with wavelength, the residuals in frames e and / will be positive. 



Fig. 6. Distribution of T and p values estimated with IRIS 
100 yum data and Planck data at frequencies 353 GHz and higher. 
The values have been calculated pixel-by-pixel using maps at a 
resolution of 5'. Left frame contains all pixels inside the cover¬ 
age of SPIRE maps. Right frame contains from each field only 
25% of the pixels with the highest column densities (fourth quar¬ 
tile). The column densities have been estimated using the same 
data and a fixed value of spectral index, yS =1.8. 


where, like in G300.86-9.00, the 217 GHz excess is positive but 
anti-correlated with column density. Field G82.65-2.00 contains 
a long filament that is visible as a region of lower p and smaller 
217 GHz excess. In the densest parts of the filament the excess 
is approximately zero. In G6.03-r36.73 (LDN 183), the high col¬ 
umn density clump is also associated with a low 217 GHz excess 
and slightly loweryS. In these fields, the temperature variations in 
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Fig. 4. T and p values for the 116 fields, calcu¬ 
lated from IRIS and Planck data (v > 353 GHz) 
at FWHM=10' resolution. The thin contours 
show two examples of the error regions, cover¬ 
ing 68% (red contours) or 95% (blue contours) 
of the distributions. The thick contours corre¬ 
spond to 50% smaller uncertainties (see text). 
The dashed c yan and magenta line s sho w T{p) 
relation from|Dupac et al.| ( |2003| ) and |Desert| 
|et al.| ( [2008] ), respectively. Frames b and c show 
the marginal distributions of p and T. Frame d 
shows the distribution as a colour image (for 
easier comparison with Fig.[^. 


the densest clumps may lower the apparent value of the spectral 
index so much that the long-wavelength excess disappears. We 
return to the question of the long wavelength excess in Sect. |4.4[ 
in particular regarding the uncertainty of the CO correction. 

In a few fields, the 217 GHz residuals show a local excess 
that may be caused by free-free emission from HII regions (e.g., 
G92.63-10.43, G107.20-f5.52). These are relatively easy to sep¬ 
arate from the extended emission where the absolute 217 GHz 
excess follows the column density structure. The extended free- 
free emission should already be accou nted for by the subtraction 
of the low frequency model (see Sect. |2.2| ). 


3.2. General results based on SPIRE data 


To gain access to smaller spatial scales, we must turn to Herschel 
data. The observations extend only to 500yum, making them less 
efficient in constraining the spectral index. We start by examin¬ 
ing the colour temperature and spectral index values determined 
with SPIRE data only. As shown in Fig. the parameters T 
and p are strongly degenerate, i.e., the colour variations caused 
by temperature and spectral index changes follow almost par¬ 
allel lines. To measure AyS ~ 0.2 variations with these bands, 
the surface brightness data must have an accuracy of almost 1%. 
Nevertheless, we analyse these data in an attempt to obtain re¬ 
sults that are mostly independent from the previous IRAS and 
Planck fits. This is particularly relevant because IRAS lOOyum 
data have been used by most studies where a negative correla¬ 
tion between colou r temperature and spectral index has been re¬ 
ported (for examplepupac et al. [2003[[Desert et al. [2008 [ [Planck 


[Collaboration XI|26l4[ ). Our only dependence on IRAS is in the 

250 yum band where the zero point and gain correction depend on 
interpolation between IRAS 100yum and Planck 857 GHz bands. 
Even there the weight of IRAS data is relatively small, thanks to 
the large distance in wavelength. 

The relative calibration of SPIRE channels is estimated to 
be better than 2%|^ However, when the gain calibration was ad¬ 
justed using Planck and IRIS data, the formal uncertainties were 
smaller. The main concern is the accuracy of the zero point cor¬ 
rection (see Sect. 2.1), which has a large impact at low column 
densities. The distributions of the estimated statistical errors of 
the intensity zero points are shown in Fig. The median val¬ 
ues are 0.24, 0.13, and 0.05 MJy sr“^ at 250yum, 350yum, and 
500 yum, respectively. For a surface brightness of 10 MJy sr“^ the 
statistical errors are thus only at the level of 1%. Frame b shows 
the distribution of relative errors calculated using the median sur¬ 
face brightness of each field. For the lowest column densities, the 


^ SPIRE Observer’s manual, 
httpiUherschel. esac. esa. int/Documentation. shtml 



Fig. 7. Illustration of the degeneracy between T and yS in a 
SPIRE colour-colour plot. The solid lines correspond to constant 
values of p or T. The lengths of the lines correspond to a tem¬ 
perature range of 10 K-23 K or a spectral index range of 1.0-2.5, 
respectively. Four distinct parameter pairs are plotted as red dots. 
The surrounding contours indicate the FWHM extent of the es¬ 
timated parameter distributions in case of 1% or 5% uncertainty 
in all surface brightness data. 


relative errors can be larger by a factor of a few. If the zero point 
errors are correlated between the SPIRE bands, the dispersion 
of the derived T and p values would be smaller. In other words, 
SPIRE data should be accurate enough to see definite trends, 
although the absolute values may be affected by systematic er¬ 
rors. Because the offset and slope of the linear fits {Herschel vs 
reference data) are anticorrelated and the zero point and calibra¬ 
tion errors are added in squares, the total estimated uncertainty 
of SPIRE measurements may be overestimated. The effect be¬ 
comes negligible outside the faintest regions {S (250yum) above 
~20MJy sr“^) where the errors are dominated by the multiplica¬ 
tive component. 

To study the general behaviour of colour temperature and 
spectral index values and to examine the effect of gain and zero 
point uncertainties, we examined pixel samples that were se¬ 
lected in each field using a uniform criterion on surface bright¬ 
ness. It was already known that the highest column densities 
are associated with the lowest colour temperatures (estimates 
based on a fixed value of P). We selected pixels where the 
350yum surface brightness was within 10% of a selected sur¬ 
face brightness level Sq. The use of 250yum surface brightness 
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Fig. 8. Distributions of estimated statistical errors cr in SPIRE 
surface brightness zero points for our fields. The values indi¬ 
cated are the formal errors of linear fits between Herschel and 
the interpolated IRIS and Planck data. The estimated error in 
the 250yum zero point is larger than 1.3 MJy sr“^ for three fields. 
Frame b shows the same distribution as a function of relative 
uncertainty, cr divided by the median surface brightness of the 
map. 


would result in a very similar selection because the ratio be¬ 
tween 250yum and 350yum surface brightness changes only by 
-13% in the temperature range 15+2 K (assuming a fixed spec¬ 
tral index). For a temperature of 15 K and {3 = 1.8, the ratio is 
/(250yum//(350yum) -1.6. In each field, we then calculated the 
average surface brightness values of the selected pixels in the 
three SPIRE bands, and used this restricted set of measurements 
to estimate T and p. The motivation for this approach is that 
when surface brightness is averaged over a large number of pix¬ 
els, the statistical uncertainty becomes very small. The errors are 
then dominated by calibration and zero point errors, and poten¬ 
tially by other systematic errors that affect entire maps (or parts 
of them). 

Figure shows the results as a function of poten¬ 
tial errors in the 250yum data. The diffuse medium with 
5'(350yum) -lOMJy sr“^ (red band) shows higher temperatures 
and lower spectral index values than the denser regions with 
5'(350yum) -lOOMJy sr“^ (blue band). The temperature differ¬ 
ence of AT = +5 K corresponds to AyS of about -0.2. The tem¬ 
perature difference of AT = +5 K corresponds to AyS of about 
-0.2. The left panels of Figure]^ show that errors in the gain cal- 



Fig.9. The T and p values for pixels with 350yum surface 
brightness -lOMJysr"^ (red band) and -lOOMJy sr“^ (blue 
band). The values are plotted as a function of multiplicative (left 
frames) and additive (right frames) modifications of the mea¬ 
sured 250yum surface brightness. The solid curves are the median 
values (calculated over all fields) and the shaded region between 
the 25% and 75% percentiles indicates the field-to-field disper¬ 
sion. The grey vertical bands correspond to the probable uncer¬ 
tainties of -2% in the relative gain calibration and ±0.2 MJy sr“^ 
in the surface brightness zero point. 



Fig. 10. The same as Fig. but for hypothetical errors in the 
350yum surface brightness. 


ibration would affect the absolute values of T and p, but would 
leave the difference between the two pixel samples almost un¬ 
changed. The right panels, by contrast, show that zero point 
errors at a level of 0.2MJysr“^ would have little effect on the 
5'(350yum) -lOOMJysr"^ sample, but a much larger impact on 
the S (350yum) -10 MJy sr“^ sample. Nevertheless, the difference 
between the two pixel samples disappears only if we assume that 
most fields are affected by a systematic error of -0.6MJysr“^ 
This would be larger than the typical estimated statistical uncer¬ 
tainties. 

The situation changes somewhat if we consider similar un¬ 
certainties in the 350yum band (Fig.[^. The multiplicative er¬ 
rors (left panel) again have little importance but the results 
are sensitive to zero point errors. A shift of A5'(350yum) = 
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Fig. 11. The same as Fig. but using surface brightness data 
after subtracting the local background. 


0.2MJysr“^ is enough to make {3 values equal in the two sam¬ 
ples. A larger offset would result in the low surface bright¬ 
ness pixels having larger p values than the high surface 
brightness pixels. If statistical errors of the zero points were 
—1-0.2 MJy sr“\ the scatter of p values should be twice as large 
as what is actually observed. For example, for the lOOMJy sr“^ 
sample alone, the apparent p values would be expected to range 
between 1.7 and 2.3 (1 cr range; see Fig. [T^. If the errors are 
statistical and without significant bias. Fig. |1()| still suggests that 
P is positively correlated with column density. However, to be 
confident of this general conclusion, systematic zero point errors 
would need to be constrained to a level of ~0.1 MJy sr“^ 

We can derive alternative T and p estimates where the zero 
point errors are completely independent of those in Fig.[^ This 
is possible by subtracting the local ba ckground. We use th e ref¬ 
erence regions defined in Table 1 of Juvela et al. ( |2012| ). The 
results are shown in Fig. The selected 350 fim surface bright¬ 
ness levels are the same as above, 10 MJy sr“^ and 100 MJy sr“^ 
Because these are now background subtracted values, they corre¬ 
spond to somewhat higher values of absolute surface brightness, 
the actual difference depending on the field. Nevertheless, the re¬ 
sults are remarkably similar to Fig. The difference between 
the lOOMJysr"^ and lOMJysr"^ samples is AT = -4K and 
Ap = -hO.2. As before, if the 350yum data contained a system¬ 
atic error of -0.2 MJy sr“\ the data would be consistent with no 
spectral index variations correlated with surface brightness (or 
temperature). However, this would require the same systematic 
error in most fields. While the zero point correction with Planck 
and IRIS data could contain some bias, a similar systematic er¬ 
ror is more difficult to understand in the case of local background 
subtraction that is carried out independently in each field. 

In summary, the SPIRE data gives some indication of spec¬ 
tral index variations, with higher p values being found towards 
regions of higher column density and lower temperature. The re¬ 
sults are made more plausible by the consistency between the 
independent analyses conducted using absolute and relative zero 
points. Because SPIRE bands cover a relatively narrow range of 
wavelength, data must have very high precision for the absolute 
values of T and p to be measured. SPIRE channels may not be 
enough for accurate determination of the absolute values of T 
and p. Nevertheless, the above experiments show that these data 


have very low systematic errors and, therefore, should provide 
reliable results when combined with data at other wavelengths. 


3.3. Genera/resL//fs from Herschel 160jam-SOOjam data 


The inclusion of PACS 160 jam data helps to constrain the colour 
temperature and should therefore also yield more accurate esti¬ 
mates of spectral index variations. This advantage is diminished 
by (1) the smaller size of the PACS maps, (2) the uncertainty of 
the relative calibration between PACS and SPIRE, (3) a possible 
bias resulting from the greater sensitivity to the warmest dust 


along the line-of-sight (|Shetty et al.|2009b 

Malinen et al.|2011l 

Juvela & Ysard|2012bl|Pagani et al.|2015 

, and (4) a larger un- 


certainty in the determination of the zero point. The last point 
results from the smaller size of the maps, and from the lower 
signal-to-noise ratio at shorter wavelengths for cold dust emis¬ 
sion. Points (3) and (4), together with the contribution of VSGs, 
are the reason why we refrain from using the 100 jam PACS data 
in our analysis. 

We continue to use a 2% uncertainty for the gain calibration 
in the SPIRE bands (band-to-band accuracy) and a 10% uncer¬ 
tainty for the 160 jam measurement^ the latter number also cov¬ 
ering the uncertainty in the relative calibration of instruments. 

In spite of its problems, the inclusion of 160yum data results 
in smaller formal uncertainties of dust temperature, especially 
when estimated using a constant value of p. Thus, in the follow¬ 
ing we compare data sets that are selected based on dust tem¬ 
perature or dust column density instead of using surface bright¬ 
ness thresholds (as in Sect. |3.2| ). We start by looking for global 
trends in p in relation to dust temperature. We selected in each 
field all pixels where the colour temperature was within 20% of 
T = IdKorT =18 K, based on the fits made with a fixed value 
of p = 1.8. Fields that did not contain at least 20 pixels in both 
categories were excluded. 

Figure shows the median values of T and p and their sen¬ 
sitivity to surface brightness errors. As seen for the SPIRE data 
alone, the p values are higher for the cold pixel sample, although 
the field-to-field scatter is comparable to this difference. Before 
the 160yum zero points were adjusted using the AKARI compar¬ 
ison, the mean p value of the warm pixels was actually above 
the average p of the cold pixels. This illustrates the sensitivity 
of our analysis to zero point errors, at least for pixels with low 
brightness (i.e. at the lOMJysr"^ level). If we use local back¬ 
ground subtraction, the situation remains essentially the same, 
with the cold pixels now showing p values that have typically 
increased by ~ 0.1. Because of the smaller size of the 160yum 
maps, we cannot use the previous reference areas. Therefore, 
the background was determined by selecting pixels with surface 
brightness between the 1% and 20% percentiles, and calculating 
the average over the subset of those pixels that is common to all 
four bands. 


Systematically higher p values are seen also in regions of 
higher optical depth. In Fig. we compare pixels within 20% 
of the values T(250yum) = 1 x 10“^ and T(250yum) = 4x 10“^. The 
optical depths were derived with a fixed value of/5 = 1.8 and are 
thus not affected by the noise-induced anti-correlation between 
T and p that could lead to some artificial correlations. Pixels in 
the higher optical depth interval have ~ 0.1 units higher p val¬ 
ues. The figure shows the results obtained with local background 
subtraction. Using the estimated absolute zero points of surface 


^ https://nhscsci.ipac.caltech.edu 
/sc/index.php/Pacs/AbsoluteCalibration 
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Multiplication of S(250//m) 


Offset for S(250//m) (MJy sr“^) 


Fig. 12. Comparison of median colour temperature and spectral 
index values derived from Herschel 160yum-500yum data for all 
pixels where modified blackbody fits with yS = 1.8 yield colour 
temperatures within 20% of 14K (solid lines) or 18K (dashed 
lines). The lines show the dependence on potential multiplicative 
or additive errors in individual bands. The red, green, blue, and 
black colours correspond to the four bands in order of increasing 
wavelength. Assuming no errors (centre of x-axis), the median (3 
is ~ 0.1 higher for the sample of cold pixels. The shaded regions 
indicate the field-to-field variation (interquartile ranges). 
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but comparing pixels with 


,-3 


Fig. 13. Similar to Fig. 

T(250/im) ~ 1 X 10“^ (dashed lines) and T(250yum) ~ 4 x 10 
(solid lines) and using T and p values calculated from surface 
brightness after the subtraction of local background. 


brightness, the results are similar, the difference of p values be¬ 
ing -0.08 instead of -0.10. 

Thus, systematic changes in p are seen for pixels selected 
based on either temperature or optical depth. Part of the system¬ 
atic effects may be masked by the field-by-field variations of the 
local radiation field. Also, the strongest changes may be limited 
in small regions within each field. Therefore, we need to inves¬ 
tigate the p variations also directly on maps. 


3.4. Herschel T and p maps 


The temperature and spectral index maps of all 116 fields are 
shown in Appen dix [C| These were calculated with a basic 
MCMC method ( |Juvda et al.||2013 ) using either 160-500yum 
data (four bands) or only data with 250-500yum (three bands). 
The error estimates are based on the assumption of 2 % and 10 % 
relative uncertainty for the SPIRE and PACS channels, respec¬ 
tively. These estimates (converted pixel-by-pixel to absolute un¬ 
certainties of the surface brightness) and the estimated absolute 
uncertainties of the surface brightness zero points were added 
together in quadrature. The PACS error estimate covers the un¬ 
certainty of the relative calibration of the two instruments. 

The appearance of the T and p maps and the dispersion be¬ 
tween neighbouring pixels gives some impression of the purely 
statistical noise but they are not useful for estimating systematic 
effects that zero point errors may cause. In the worst case, a zero 
point error can even change the sign of the T-p correlation. We 
calculated 100 realisations, varying the surface brightness zero 
points within their estimated uncertainties, each calculation re¬ 
sulting in different T and p maps. The realisations were sorted 
according to the average values of p. The p maps correspond to 
the 16% and 84% percentiles of the distribution and 
were plotted and can be used to indicate how sensitive the mor¬ 
phology of the T and p maps is to zero point errors. The p^^"^ 
and p^^^ maps are estimated through minimisation. Because 
of the a priori constraint on T and p, the level of the correspond¬ 
ing MCMC p maps can be slightly different. 

The maps derived from SPIRE data are shown in Appendixj^ 
and two examples are included below (Figs.[^and[^. An anti¬ 
correlation between T and p is seen in most fields. There are 
some twenty fields for which the p variation is not clearly visible 
above random fluctuations or map artefacts. In a few cases (e.g. 
G110.89-2.78 and G163.82-8.44), the observed large-scale anti¬ 
correlation is no longer significant when the zero point uncer¬ 
tainties are taken into account. Nevertheless, an anti-correlation 
is detected in nearly 90% of the fields, showing the ubiquity of 
this behaviour. However, in a given field, the effect can be re¬ 
stricted to a small area or even a single core. 

There are a few fields where the data show some positive cor¬ 
relations between T and p, mainly in small regions near strong 
point sources (e.g. in fields G343.64-2.31 and G345.39-3.97). In 
G82.65-2.00, which is shown in Fig.[^ the correlation between 
T and p is also partly positive, this being associated with a cold 
filament that runs diagonally across the field. At smaller scales, 
at the level of individual clumps, the parameters are still pref¬ 
erentially anti-correlated. In Sect. |4.5[ we argue that a positive 
correlation could in some cases be explained by radiative trans¬ 
fer effects that occur when the clouds have high optical depth for 
the radiation heating the dust. 

When PACS 160yum data are included, the results on the T- 
P correlations are qualitatively similar although the values of¬ 
ten differ. The four-wavelength fits have lower formal error es¬ 
timates but are affected by the larger uncertainty of the relative 
calibration of the two instruments and the different methods used 
in the zero point correction of PACS and SPIRE data. The me¬ 
dian values over all fields, calculated over PACS coverage, are 
T = 16.9KandyS= 1.89 for the 250-500yum fits and T = 18.4K 
and p = 1.66 for the 160-500yum fits. In the latter case, part of 
the higher temperature may be attributed to the fact that shorter 
wavelengths are more sensitive to line-of-sight temperature vari¬ 
ations and especially to the warmest regions. However, consider¬ 
ing the other sources of uncertainty, this cannot be concluded to 
be necessarily the main source of the differences. Based on the 
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Fig. 14. An example of negative correlation between colour tem¬ 
perature and spectral index. The field is G300.86-9.00 in Musca. 
The frames, as indicated by their colour scales, show 250yum 
surface brightness, t(250 /im) optical depth, colour temperature 
r, and the spectral index (3. The bottom frames show the p maps 
that correspond to the 16% and 84% percentiles in a Monte Carlo 
study of zero point errors (see text). The maps are based on 
SPIRE data and, for plotting only, they have been convolved to 
a resolution of 1'. 
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Fig. 15. Field G82.65-2.00 as an example of some large scale 
positive correlation between colour temperature and spectral in¬ 
dex. The frames are as in Fig.p^ 


models shown in Appendix [A| the inclusion of the 160yum band 
could be expected to lead to changes AT -0.2 K and AyS —0.04 
(in case of the ry = 10 mag model and compared to 250-500 yum 
data). 


3.5. Combined analysis of Herschel and Planck data 

A final set of T and p maps was calculated using Herschel 250- 
500yum data and the four Planck bands between 857 GHz to 
217 GHz. The goal is to derive the most reliable estimates of 
P, using the high resolution Herschel observations while at the 
same time making use of Planck data to constrain p at large 
scale s. Th e Herschel data are kept at the ~ 37" resolution (see 
Sect. \2.\) . The 160yum maps are not used, mainly because of 
their smaller size. The resulting maps of T and p are shown in 
Appendix [D| 

In these calculations the observed data could be used at their 
original resolution but we continued to use Planck maps con¬ 
volved to a 5.0' resolution. The 217 GHz data are used at the 
original 5.01' resolution, the pixel sampling being visible in the 
result ma ps. T he Planck and IRIS data were colour corrected 
(see Sect. 133 but the colour corrections were not updated dur¬ 
ing the fitting because their effect on surface brightness values is 
very small. Furthermore, if the spectral index is not constant over 
the full wavelength range, it is not clear that the joint fit would 
result in a better description of the SED shape over a given band. 

The model is defined on 28" pixels and it consists of one 
intensity, colour temperature, and spectral index value per pixel. 
The;^^ values that are minimised consist of two components. For 
Herschel bands, the model is compared to the 250-500yum data 
directly, treating each pixel as an independent measurement. The 
second component of the value results from the comparison 
with Planck data, requiring the convolution of the model predic¬ 
tions to the 5.0' resolution. This makes the calculations computa¬ 
tionally demanding, partly because of the convolution operation 
itself but mainly because of the correlations that the convolu¬ 
tion operation introduces between model pixels. Nevertheless, 
the procedure enables us to seek a solution at Herschel resolu¬ 
tion that is still constrained by the Planck data on larger scales. 

The total probability of a model is a combination of the 
probabilities of fits to three Herschel and four Planck bands. 
Considering that each Herschel map has been rescaled using 
the same Planck and IRIS data, we adopt for SPIRE channels 
an error estimate of 3.5%. This is a compromise between the 
original 2% band-to-band uncertainty and the 7% uncertainty of 
the absolute calibration. The zero point uncertainties are added 
to these numbers in quadrature (see Sect. |3.4[ ). For the Planck 
data, we keep the 10% error estimates for the 857 GHz and 
545 GHz channels and 2% at lower frequencies. The 353 GHz 
and 217 GHz error estimates were further increased with an 
amount corresponding to 50% of the CO correction. Because of 
this, the relative error estimates are often largest in the 217 GHz 
band. 

Because Herschel has many independent measurements over 
the area covered by a single Planck beam, weighting according 
to the above listed error estimates gives Herschel a large weight 
compared to Planck. However, at large scales the systematic er¬ 
rors should be smaller for Planck because the data are taken from 
continuous all-sky maps. Herschel maps can suffer from minor 
gradients (e.g. due to map-making) that cannot be recognised be¬ 
cause of the small size of the maps. Therefore, we increased the 
relative weight of Planck data in the;^^ calculation by an ad hoc 
factor of three. This is not a critical issue because the actual fit 
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Fig. 16. Modified blackbody fits in field G300.86-9.00 using the combination of Herschel and Planck data. The uppermost frames 
show the fitted intensity at 250yum, and the colour temperature and spectral index maps. The relative residuals (observation minus 
model, divided by model prediction) are shown in frames d and / for the Planck bands of 353 GHz and 217 GHz. Frame / shows 
the median residuals for the three Herschel bands (red symbols) and the four Planck bands (blue symbols). The lower dotted line 
corresponds to twice the default CO correction. The error bars correspond to the median error estimate of the surface brightness 
over the map. 


only provides a convenient point of reference, the detailed infor¬ 
mation being contained in the fit residuals. 

The calculated T and {3 maps are shown in Appendix [P] One 
example is shown in Fig. 16 The T and (3 maps (frames b and 
c) are relatively similar to the earlier Herschel fits, again show¬ 
ing clear anti-correlation between the parameters. The area in 
Fig.[^is ~80 times the solid angle of a 5' Gaussian. Thus, con¬ 
sidering the potential effects of noise, the maps still contain a 
large number of independent samples. Frame / shows the me¬ 
dian residuals for the seven bands used in the fits. The residuals 
are small in all Herschel bands, indicating that these still carry 
most weight in the overall fit. The fit is also relatively good in 
the Planck 857 GHz-545 GHz bands, the wavelength range that 
is common with SPIRE. Although Planck shows systematically 
slightly higher surface brightness values, the residuals are at a 
level of one per cent and do not show a strong pattern that would 
be directly related to the column density structure. However, in 
the 353 GHz band (850yum) and especially in the 217 GHz band 
(1380/im) the residuals are clearly positive (see frames d and e). 
In this case, a single modified blackbody that fits 250 jam-550 jam 
data is not a good description of millimetre emission. The lower 
curve corresponds to the larger CO correction that is twice the 
default correction (i.e., assumes line ratios (2-l)/(l-0)=1.0 and 
(3-2)/(l-0)=0.6, in units of Kkm/s). In the field G300.86-9.00, 
this would be sufficient to make residuals consistent with zero. 
However, over the area covered by SPIRE maps, Planck Type 2 
CO maps indicate a line ratio of -0.4 (a least squares fit gives 
a slope of 0.39 and the direct ratio of the two maps gives 0.44 
for the average of all pixels and 0.41 for the brightest pixels in 
the fourth quartile of IT(1 - 0)). The ratio is thus lower than the 
assumed default ratio of 0.5. This suggests that the true residuals 


Table 2. The median values and 1 cr scatter of T and [3 calculated 
for all fields. In each field, the estimates correspond to spatial 
averages over a FWHM=10' beam or, in the case of first row, 
FWHM=30'. The values in parentheses indicate the 1 cr scatter 
between the fields. 


Data 

r(K) 

p 

IRIS + Planck, FWHM=30' 

16.92(1.36) 

1.71(0.10) 

IRIS + Planck, no 217 GHz 

16.03(1.31) 

1.84(0.11) 

IRIS + Planck, with 217GHz 

16.43(1.59) 

1.75(0.13) 

SPIRE 

14.89(1.25) 

2.04(0.16) 

SPIRE + l6Qpm 

17.46(2.09) 

1.68(0.15) 

SPIRE + Planck, with 217 GHz 

14.93(1.16) 

2.03(0.12) 


might be even slightly larger than indicated by the upper curve 
of Fig.p^. 

4. Discussion 

4.1. The main results and their robustness 

Table [2] summarises some of the main conclusions of our anal¬ 
ysis: (i) (3 values within the fields are larger than the average 
value at larger scales (see Fig. |^, (ii) the (3 values are smaller 
when the fits are extended to 217 GHz, and (iii) the inclusion of 
PACS 160yum data results in higher [3 values than an analysis 
using the SPIRE bands only. The first two facts are likely to re¬ 
flect real physical changes, such that high density regions exhibit 
with steeper emission spectra. On the other hand, the difference 
between the SPIRE and SPIRE-r 160 pm results may indicate the 
presence of small systematic errors. 
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Figure [^showed that there is a clear difference between the 
10' and 30' scales. Because the larger beam also covers diffuse 
regions around our target clouds, this already suggests that p is 
correlated with column density. More direct evidence of depen¬ 
dence between dust temperature, column density, and spectral 
index is provided by the maps presented in Appendices [B]|D| In 
Fig.j^yS values were ~ 0.1 units lower when the 217 GHz band 
was included in the fits, the exact value depending on the relative 
weight given to the different bands. In Sect. |3.5| and Appendix[D| 
the 217 GHz residuals are seen to be clearly positive with only 
a few exceptions. The wavelength dependence of p is examined 
further in Sect. 14.41 


clearly also in diffuse regions where the relative importance of 
the CO correction should be smaller. 


4.2. T andp correlations between fields 


In Fig. p7p9| we examine the colour temperature and spectral 
index values calculated using the average surface brightness of 
each field. Fig. [T^ shows the (T, P) relations derived with (1) 


The T and p maps were estimated with Planck data 
(Appendix [B]), with Herschel data (Appendix [C]) data, and with 
joint fits of both Planck and Herschel data (Appendix [D|). The 
anti-correlation between T and p cannot be the result of pure sta¬ 
tistical noise because that would not result in spatially correlated 
variations in the maps. The only possible explanation would be 
some kind of a systematic error that affects both Planck and 
Herschel data in a similar fashion. This is technically possible 
because the SPIRE data were scaled and zero-point corrected 
using IRIS and Planck data. However, in Sect. we showed 
that zero point errors are important only at low column densi¬ 
ties. Furthermore, the SPIRE-only results remained essentially 
the same when we used local background subtraction instead of 
absolute zero points. We also demonstrated that gain calibration 
errors are not able to significantly affect \hoT - p rel ation , not 
even when the data consists of SPIRE bands only (Sect. |3.2[ ). The 
Herschel results were relatively similar to our previous analysis 
where the SPIRE gain calibration was not tied to Planck. The fi¬ 
nal possible error source are mapp ing errors that cou ld produce 
artificial T -p anti-correlation (e.g. Juvela et al.|2QI3| ). Mapping 
artefacts are difficult to detect in surface brightness data even 
when their effect on T and p estimates should be easily recognis¬ 
able. The maps in Appendix are thus a good indicator of this 
kind of errors. Artefacts would appear as spatial features that 
are different between the 160-500/im and 250-500yum fits (es¬ 
pecially regarding artefacts in PACS data) and/or uncorrelated 
with surface brightness. As an example, one can look at the data 
on the field LDN 134 (G4.18-r35.79) in Appendix |C There are 
clear features that are following changes in column density and 
are similar in the two band combinations. However, the P/^ map 
also exhibits a north-south gradient that is correlated with neither 
the column density nor the morphology of the P'i map. This sug¬ 
gests that the PACS map suffers from a surface brightness gra¬ 
dient (possibly caused by the strong emission that extends over 
the southern map boundary and thus complicates map-making). 
Smaller scale artefacts are clear, for example, in the south corner 
of the Ts andmaps of the field G126.63-r23.55. However, in 
most cases the two band combinations give comparable results, 
especially for compact structures, and the observed p changes 
are well correlated with column density variations. This remains 
true even in the case of both G4.18-h35.79 and G126.63-r23.55. 
We therefore conclude that T - p anti-correlation is a common 
physical phenomenon, which we robustly detect in most of our 
fields. 

Regarding the change of p at wavelengths longer than 
500 yum, the result could be affected mainly by errors in the 
correction for CO emission. However, even in the joint fits of 
Herschel and Planck data, the 217 GHz excess is, after the re¬ 
moval of the expected CO contamination, typically more than 
20% of the total signal. This is large compared to the expected 
CO contamination and is thus very unlikely to be caused by inac¬ 
curacy of the CO correction. Eurthermore, the excess is detected 


IRIS and Planck, (2) SPIRE, and (3) SPIRE and Planck. The 
217 GHz Planck band was omitted from these fits. The error bars 
include statistical and calibration errors, including at 353 GHz 
the error associated with the CO correction. Eor the SPIRE chan¬ 
nels we used an error estimate of 7%. All three frames show 
some degree of T-p anti-correlation but there are also differ¬ 
ences. Eor example, the joint fit of Planck and Herschel data 
results in a wider temperature distribution. When Herschel data 
are used, the error bars strongly overestimate the scatter. This is 
particularly clear in the middle frame. Thus, Eig.p^ shows that 
the procedure for the re-scaling of Herschel data (including the 
zero points) has been fairly accurate. This also partly justifies 
the lower 3.5% relative uncertainty used for SPIRE channels in 
Sect.p31 

Figure shows similar plots using pixels in given temper- 
ature intervals. The pixel selection is based on T and T(250yum) 
values calculated with SPIRE data and a constant value of 
P = 1.8. Therefore, the selection is not directly affected by 
the T-p anti-correlation (whether physical or noise-produced). 
In Eig. the three temperature intervals are clearly separated 
but each forming a narrow distribution that extends over a wide 
range of p values. The overall dispersion might be caused by 
noise but also by the fact that fory(5 = 1.8 a given colour temper¬ 
ature interval actually contains different p values. When p was 
allowed to vary freely, some points have moved to higher tem¬ 
perature, especially from the lowest ~ 14 K bin. The fact that the 
orientation of these elongated clouds of points is not the same 
as that of the general T-p anti-correlation does not directly indi¬ 
cate which of the two effects is stronger, i.e., to what extent the 
anti-correlation is caused by noise. However, the comparison of 
the median values of the three temperature samples shows the 
anti-correlation, even when the samples were selected based on 
temperatures calculated with a fixed value yS = 1.8. 

Eigure[^ compares values calculated for three column den¬ 
sity intervals. The (T, P) distributions of the three optical depth 
intervals overlap but the median p increases with column density. 
The Pn^.T slope also becomes flatter with increasing T(250yum) 
for all three combinations of surface brightness data. Because 
the x-axis values have been calculated with a fixed value of p, 
the noise-induced anti-correlation should be reduced. However, 
higher T(250yum) values usually correspond to higher signal-to- 
noise ratios and the flattening could be partly a noise effect. In 
the case of SPIRE data (middle frame), all relations are steeper, 
this also being possibly a noise effect. The p values are not quite 
systematically higher for the higher column densities, although 
the lowest column density sample does correspond to the high¬ 
est temperatures and the lowest values of p. While Eig. paused 
the average surface brightness of all pixels within a given op¬ 
tical depth range, Eig. 19 examines each held separately. This 
increases the risk that result become affected by local artefacts 
(e.g, near map boundaries). Depending on the range of column 
densities within a map, each plotted symbol may correspond to 
a very different number of individual map pixels. 

Spatial averaging of surface brightness data should increase 
the bias caused by the mixing of spectra with different colour 
temperatures. Eigure compares T and p estimates derived 
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Fig. 17. The (r,/5) values based on the average 
surface brightness of the fields. The values are 
based on Planck 857 GHz-353 GHz and IRIS 
100 jjm data (left frame), SPIRE channels (cen¬ 
tre frame), or all the previous bands combined 
(right frame). 
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Fig. 18. As in Fig. [^but using the average surface brightness 
in three temperature intervals that are determined with SPIRE 
data and a fixed value of the spectral index, (3 = 1.8. The blue, 
green, and red symbols correspond to (14+0.5) K, (16+0.5) K, 
and (18+0.5) K, respectively. Each point corresponds to one 
field. The black lines connect the median values of the three tem¬ 
perature samples. 
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Fig. 19. As previous figures but using average surface bright¬ 
ness values in three T(250/im) intervals. The blue, green, and 
red symbols correspond to T(250yum) equal to (1 + 0.5) x 10“^, 
(2 + 0.5) X 10“^, and (3 + 0.5) x 10“^ respectively. The optical 
depths are derived from a fit to the SPIRE data using a fixed 
value of p = 1.8. The lines are least-squares fits to the points of 
the same colour and the slope values k are listed in the frames 
(unit K“^). The median values are indicated with crosses (in the 
middle frame the red and green crosses overlap) and the median 
values of p are also given. 



T(K) 


Fig. 20. Effect of spatial resolution on field-averaged parame¬ 
ters. The grayscale image shows the distributions of {T,p) values 
calculated from IRIS and Planck data for pixels. The contours 
show the distribution separately for pixels below (blue contours) 
and above (red contours) the average column density in each 
field. The contours are drawn at levels of 20%, 50%, and 90% 
of the maximum value (data points per area) for distributions 
calculated from surface brightness data at 5' (solid contours) or 
10' (dashed contours) resolution. The dotted contours (50% level 
only) correspond to a case where the surface brightness values 
of all pixels below or above the median value (per field) are av¬ 
eraged before {T,p) calculation. 


is no systematic shift in p values. This suggests that the effects 
of temperature mixing are not very pronounced at scales above 
5', up to scales approaching the size of individual fields. 

4.3. Correlation with Galactic location 


from IRIS and Planck surface brightness maps at a resolution 
of 5' and 10'. The data were further divided to low and high 
column density parts using column densities calculated using a 
fixed value yS = 1.8. As a third alternative, we averaged the sur¬ 
face brightness of all pixels above and below the median column 
density, thus calculating two T and p values per field. This last 
alternative is expected to maximise the bias and should lead to 
smaller values of the spectral index. The two column density 
samples are displaced relative to each other, the high column 
density pixels extending towards lower temperatures and higher 
values of the spectral index. However, the differences between 
the 5' and 10' resolution data are very small. In particular, there 


Figure!^ shows {T,p) values for fields in different distance in¬ 
tervals. The fields with low T and high p values are more likely 
to be located within 500 pc. For individual clumps a dependence 


on distance could be explained by beam dilution but Fig. 21 is 
concerned with the average surface brightness of entire fields. 
Nevertheless, in our sample the nearby fields are more likely to 
contain a single cold clump that dominates the field-averaged pa¬ 
rameters. Thus, the dependence between p and distance is likely 
to be mainly a selection effect. 

Figure shows the correlation with Galactic latitude. The 
scatter is again large but high latitude clouds are preferentially 
associated with higher p values. The origin of this correlation 
is partly the same as in Fig. high latitude clouds also be¬ 
ing nearby. However, compared to the case of different optical 
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Fig. 21. As previous figures but with colours indicating the esti¬ 
mated distance d of the field. The blue, green, and red symbols 
correspond to distances d < 200 pc, 200 pc < d < 500 pc, and 
d > 500 pc, respectively. The median values are indicated with 
the crosses of the same colour. The fields without reliable dis¬ 
tance estimates are marked with black symbols. 



T (K) T (K) T (K) 

Fig. 22. As previous figures but with colours indicating the 
Galactic latitude. The blue, green, and red symbols correspond 
to \b\ < 5°, 5° < \b\ < 20°, and \b\ > 20°, respectively. The 
crosses indicate the corresponding median values. 


depth intervals (Fig.[T^, the median temperature is less latitude- 
dependent. 

In Fig. the dust parameters are plotted against Galactic 
longitude and latitude, the figure showing a general drop of {3 
towards low latitudes. In addition to values averaged over each 
field, T and p are indicated for pixels T(250/im) = (2.0 + 0.5) x 
10“^. Above the lowest latitudes, these are above the average 
column density and have higher values of p. The median values 
are 1.81 and 1.88 for all the pixels and for the T(250yu) ~ 2x 10“^ 
pixels, respectively. The difference disappears below \b\ ~ 4°, 
because the average optical depth rises to and above 2 x 10“^. 
The latitude dependence can indicate a change in dust properties. 
However, at low latitudes local star-formation and the long lines- 
of-sight may increase temperature variations within the beam, 
thus decreasing the apparent values of p. 

There is no clear dependence on Galactic longitude, although 
the values for / = 0-10° are marginally above the average. 
We do not see the strong syst ematic increase of P towards th e 
inner Galaxy that was found in [Planck Collaboration XI| ( [20141 ). 
Partly this is because that effect is strongest in the Galactic plane, 
whereas all our fields are at latitudes \b\ > 1° and, furthermore, 
within |/| < 30° all are above \b\ ~ 3°. Secondly, most of our 
fields are dominated by emission from within one kiloparsec, 
too near to probe physical variations related to, e.g., the Galactic 
molecular ring. 

4.4. Wavelength dependence ofp 

In the fits performed with the combination of Herschel and 
Planck data, the 353 GHz and 217 GHz data consistently show 
an excess relative to the modified blackbody that describes the 
emission at wavelengths below -500yum. In Fig. the relative 



GLON (degrees) GLAT (degrees) 

Fig. 23. Spectral index values as function of Galactic longitude 
and latitude. Values are calculated using Planck and IRIS data. 
The blue symbols correspond to the average surface brightness 
in the entire fields. Red symbols correspond to the subset of pix- 
els with T(250jum) = (2.0 ± 0.5) X 10”^ 



Fig. 24. Colour temperature and spectral index fits for surface 
brightness data averaged over FWHM=10' beams centred on 
each of the 116 fields. The data points correspond to values ob¬ 
tained from an analysis of the IRIS and Planck data (frame a), 
and of the SPIRE and Planck data (frame b). The cyan symbols 
represent single component fits that are restricted to frequen¬ 
cies y > 353 GHz. The other plot symbols indicate two com¬ 
ponent fits that employ different spectral indices below 700yum 
(red symbols) and above 700yum (blue symbols). The Planck 
217 GHz band is included in these two component fits. The plus 
signs indicate median values. 

excess was larger in diffuse than in the dense medium. This re¬ 
mains true for the joint fit shown in Fig. As mentioned in 
Sect. O this could be related to line-of-sight temperature vari¬ 
ations that decrease the apparent FIR spectral index of dense re¬ 
gions. However, the quantification of these effects would require 
detailed modelling of the individual sources. 

We examined the wavelength dependence of p using aver¬ 
age surface brightness values towards the centre of each field. 
As in Fig. the mean surface brightness is calculated with a 
Gaussian beam with FWHM equal to 10'. Figure shows the 
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single component fits for two data combinations. These consist 
of Planck data with y > 353 GHz and with either IRIS lOOyum 
band or the three SPIRE bands. For the latter combination, the 
scatter is larger but the median value of p is in both cases 1.8- 
1.9. Alternative fits were performed using two p values, one be¬ 
low and one above To = 700 /im. This model thus consists of one 
intensity parameter, one colour temperature value, and two spec¬ 
tral index values that are here named ySpiR and y^rnm- The model is 
required to be continuous at 700yum. The Planck 217 GHz band 
is included in these fits. In Fig.[^ the two component fits show 
a clear difference between the FIR part (T < 700yum) and the 
millimetre part (217 GHz-353 GHz, ~850-1400yum). In the mil¬ 
limetre part the values are y^mm = 1.6-1.7 for both data combi¬ 
nations. At shorter wavelengths, the median values are signifi¬ 
cantly higher, ySpiR = 1-96 for lRlS-\-Planck and ySpiR = 2.11 for 
SFlRE-\-Planck. 

Figure shows the SEDs that are estimated from surface 
brightness values averaged over a 10' beam towards each of the 
fields. The values are further averaged either over all fields, over 
fields with column density below or above the median column 
density, or over fields with colour temperature below or above 
the median colour temperature. The selection of the fields is 
based on values (T(250yum) and T) derived from the analysis of 
IRIS and Planck data (v > 353 GHz) but all bands are included in 
the fits shown in Fig.[^ As above, the fits employ a single dust 
temperature but two spectral index values. The derived param¬ 
eter values are given in the figure. For signal averaged over all 
fields, the spectral index values are y(5piR = 1.91 and y^mm = 1-66. 
Compared to a fit with a single p value (17.0 K and p = 1.72), 
the reduced value is smaller by a factor of three. When fields 
are divided according to column density, higher column den¬ 
sity fields are associated with slightly lower dust temperature, 
AT = 0.4 K. The values of ySpiR are practically identical for the 
two column density intervals but ySmm is 0.07 units lower for the 
lower column density sample. 

When the fields are divided based on their colour tempera¬ 
ture, differences in the spectral index are more noticeable. The 
average signal of cold fields has ySpiR that is higher by ~0.1. The 
data are averaged over 10' beams and, for both samples, over 
almost 60 fields. Thus the noise should be very small although 
the results can be affected by systematic errors if these cause 
scatter between the T and p estimates of individual fields. The 
zero point errors are one such possible error source although, by 
concentrating on high column density central areas, their effect 
should be limited. Compared to ySpiR, the parameter y^mm should 
be more insensitive to T errors but it also shows a clear differ¬ 
ence, y^rnm being 0.12 unit higher for the cold field sample. This 
is not necessarily very significant, considering the limited wave¬ 
length range. On the other hand, the 217 GHz point is in every 
case several cr above the spectrum extrapolated with ySpiR. 

The CO correction is one of the main uncertainties regard¬ 
ing the magnitude of the 217 GHz excess. Our default correction 
assumed a line ratio (2-l)/(l-0) = 0.5, which is consistent with 
the average line ratio obtained from Planck Type 1 CO maps (see 
Fig.|g. In reality the line ratio is not spatially constant and could 
be higher in individual clumps. To examine the possible conse¬ 
quences, we repeated part of the analysis using CO corrections 
that were twice as large as the default one. This corresponds to 
line ratios rA(2-l)/rA(l-0) = 1.0 and Tp,0-2)ITp,{\-Q) = 0.6. 
According to Fig. this is a very conservative estimate of the 
average (2-l)/(l-0) ratio and also exceeds the largest observed 
values of the (3-2)/(2-l) ratio. 

If the CO correction is doubled, the 217 GHz residuals of 
Fig. [^remain positive relative to the extrapolated FIR compo¬ 


nent with y^FiR- The excess is ~1 sigma for the whole sample, 
more for the low column density sample and close to zero for the 
high column density sample. In Fig. the comparison also is 
affected by spatial averaging that (via temperature mixing) tends 
to decrease t he va lue of ySpiR- We therefore also repeated the anal¬ 
ysis of Sect. |3.5| (spatial resolution ~ 5' at 217 GHz) using the 
larger CO correction. Figure [26| shows the resulting T and p es¬ 
timates and the 217 GHz excess (observed surface brightness in 
relation to the fitted modified blackbody with a single-yS value). 
The values of p again are higher in regions of higher column 
density. The difference is larger in case of the larger CO correc¬ 
tion. There also is a clear effect on T and p values. The larger 
correction leads to almost 0.1 units higher median value of p. 
The CO correction directly reduces the 217 GHz surface bright¬ 
ness attributed to dust emission and thus leads to smaller relative 
217 GHz excess. The excess is smaller within the high column 
sample, in spite of the higher p values in these regions. The ob¬ 
served surface brightness exceeds the fitted modified blackbody 
curve in ~75% of the fields. Taking into account that the larger 
CO correction very likely overestimates the average CO contam¬ 
ination, Fig. [26| strongly suggests that the 217 GHz excess is real 

and^FiR^^ mm 

Although the apparent values of ySpiR are higher, they are 
expected to be decreased by line-of-sight variations, which are 
significant in such high column density regions (see Sect. |4.5| ). 
This suggests that the intrinsic dust properties may have an 
even stronger wavelength dependence. Also, in the two compo¬ 
nent fit the exact values of p depend on the pivot wavelength. 
If this is decreased, the results remain qualitatively similar but 
values of ySpiR decrease and values of y^mm increase. For the 
To = 700yum^490yum the difference between ySpiR and yS^m de¬ 
creases by more than half. Detailed characterisation of the wave¬ 
length dependence, including the behaviour at smaller spatial 
scales, is left for future papers. 


4.5. Effects of temperature variations 


The observed p depends on the range of dust temperatures within 
the beam. The observed value will be below the intrinsic opacity 
spectral index of the dust grains, the difference increasing with 
increasing t emperature dispersion ( [Malinen et al.||20lT[ |Ysard 


et al.|[MT2 ). In our sample, we have objects with two types 


of temperature gradients: externally heated clumps and cores 
heated internally by young stellar objects (YSOs). We examined 
th e effect of YSOs using the catalog ue of submillimetre clumps 
in Montillaud et al. ( 2015, in press] ). The clumps are classified 
based on WISE photometric data and submillimetre dust tem¬ 
perature either as protostellar (with YSO candidates) or starless. 
For a large fraction of the clumps the classification remains un¬ 
determined (see [Montillaud et al.|2015, in press[ for details). In 
the area covered by both PACS and SPIRE instruments, there are 
204 protostar candidates, 1156 starless clumps, and 1968 unclas¬ 
sified sources. The distributions of the corresponding (T, P) val¬ 
ues of the clumps are shown in Fig. [27] where the spectral indices 
are those derived from 160-500yum data at a spatial resolution 
of 60". 


There is significant overlap between the three clump cate¬ 
gories and the modes of the distributions do actually coincide. 
However, starless clumps are found preferentially at lower T and 
higher p values while the distribution of protostellar clumps ex¬ 
tends to higher temperatures. The median values of p are 1.82 
and 1.71 for the starless and the protostellar clumps, respec¬ 
tively. If the YSOs heat a volume large enough to be significant 
at long wavelengths (and at the 60" resolution) the difference 
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Fig. 25. Spectral energy distributions based on surface brightness averaged over fields and over 10' beams towards the centre of each 
field. Frame a corresponds to the average over all fields, frame h to averages over fields below (blue) and above (red) the median 
column density, and frame c to averages over fields below and above the median temperature. The dotted lines are long-wavelength 
extrapolations using ySpiR- At 217 GHz the lower points (without error bars) correspond to data that have been corrected assuming a 
line ratio of (2-l)/(l-0)=1.0 and thus with twice the default CO correction. 


could be an indication either of dust processing (e.g., evapora¬ 
tion of ice mantles or dynamical processing caused by outflows) 
or of a dire ct temperature effect on the grain material ( |Meny| 
et al. |2007 ). However, the observed differences may partly re¬ 


flect the increased temperature variations that are caused by the 
YSOs, without any difference in the intrinsic p of the dust grains 


(see|Shetty et al.|2009bt|Malinen et al.|2011[[Ju^vela et al.|2011[ 

[Juvela & Ysard|2Q12b|). 

If the temperature dispersion was caused mainly by vary¬ 
ing degrees of internal heating, it could result in the kind of 
T -p anti-correlation as seen in Fig.[^ However this is unlikely 
because the effect of the (mostly weak) protostellar sources is 
highly localised and, furthermore, the anti-correlation is more 
pronounced for the starless clumps. A complementary ques¬ 
tion is whether p depends on the grain size distribution. Grains 
are known to grow inside dark clouds and their growth is re¬ 
vealed via various phenomena, most notably the coreshine ef¬ 


fect (Steinacker et al 


agani et al.|2Q10t[Steinacker et al.| 


|2014t [Lefevre et al.||2014| ). A variation of p due to the grain 
growth would result in a segregation of the coreshine cases in 
the T -yS plot (such as Figs. and [^. Paladini et al. (in prep.) 
show that this is not the case and no preferential region appears 
in their T - p relation. This gives a clue that the p value might 
not be related to the size of the grains, or more probably that 
variations in the dust composition and size, as well as tempera¬ 
ture changes along the line-of-sight do not allow us to see clear 
trends. Moreover, for starless cores with the highest values of p, 
coreshine should trace grain growth deep inside the cores while 
the r - yS relationship is not representative of the coldest region 
( jPagani et al.||2015] ). For embedded sources, only a full multi¬ 
wavelength 3D modelling, including grain property and temper¬ 
ature variations, can help to remove the degeneracy. 

In the case of externally heated clumps, the temperature vari¬ 
ations again result in the observed spectral index values being 
lower than the intrinsic opacity spectral index of the dust grains. 
However, unlike in the case of internal sources, the temperature 
variati ons are likely to favo ur a positive correlation between T 
and p ( jMalinen et al.||2QIl| ). The precise effect depends on the 
strength of the temperature gradients ( [Juvela & Ysard|2QI2b| ). 

We show in Appendix [^results from radiative transfer mod¬ 
elling of optically thick filaments. These suggest that once the 


visual extinction exceeds Ay ~ 10 mag, the observed p can be 
more than AyS = 0.1 below the intrinsic value of the dust grains. 
The effect is well correlated with the drop of temperature. When 
the opacity is increased close to Ay ~ 50 mag, the difference 
between the observed and intrinsic p values could be AyS ~ 0.7. 
However, such large decreases in p values are not seen in our 
sample of fields. For example, the main clump of G206.33-25.94 
has an estimated optical depth of T(250yum) = 0.0047 (derived 
with yS = 1.8 and 160-500yum data at V resolution). Assuming 
the average NIR-sub millimetre opa c ity rat io of T(250yum)/Tj = 
1.6 X 10“^ derived in | Juvela et al. (2015), this corresponds to 
Tj ~ 2.9 and roughly 10 mag of visual extinction. The clump is 
not associated with any dip in the p map and the absolute values 
of p are high, at leastyS ~ 2.0 (depending on the data set used). In 
field G150.47-r3.93 the peak optical depth is T(250yum) ~ 0.007, 
which corresponds to Ay ~ 16 mag. Values of p are still increas¬ 
ing towards the column density peak. 

Some drop in p is observed only towards the centre of 
G6.03-r36.73 (LDN 183) where the estimated optical depth is 
T(250yum) = 0.017 (foryS = 1.8 and a resolution of 37") suggest¬ 
ing Av ~ 50 mag. In this source, the estimated subm m dust emis- 
sivity is higher than average (see [Juvela et al.[20T5 ), but for such 
high opacities the column density derived from dust emission is 
still likely to be underestimated. Based partly on MIR absorp¬ 
tion, the estimated maximum extinction exceeds Ay = 100 mag. 


although on ly in an area with less than V in size (Pagani et al. 
2003[[2015| . At the resolution of 60", the observed p values drop 
from p ~ 2.0 in the surrounding regions down io p = 1.60 in the 
central core, the values being measured using 160-500yum data 
at a resolution of 1'. The drop is much less than would be ex¬ 
pected based on the simple models of Appendix [A| 

In Appendix we also tested cases where the dust in the 
central part of the filament is associated with higher submillime¬ 
tre opacity. This leads to lower central temperature and thus to 
larger temperature gradients. As a result, the observed value of p 
is further decreased towards the centre of the model cloud. This 
can happen even when the intrinsic p is actually higher in the 
centre. 

Compared to the models, the observations do not show such 
a strong drop of p. This could be understood if the clouds had a 
very clumpy structure (smaller temperature gradients), the high- 
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Fig. 26. Comparison of fits of Planck and Herschel data in case 
of different corrections for CO emission. The boxplots show the 
distribution of field-averaged values of T, /5, and 217 GHz ex¬ 
cess (relative to fitted modified grey body model, ~ 5' resolu¬ 
tion). Distributions are plotted dividing each field to pixels be¬ 
low (‘W < med”) and above (‘W > med”) the median column 
density. The plots on white background correspond to the default 
CO correction and the plots on grey background to a CO correc¬ 
tion that is twice as high. Each boxplot shows the median pa¬ 
rameter value (red horizontal line), the interquartile range (box), 
and the 10%-90% extent of the distribution (“whiskers”). Each 
individual field outside the [10%, 90%] range is plotted with a 
plus sign. 

est column density peaks remained unresolved, or the increase of 
the intrinsic {3 values is very large and/or extends clearly outside 
the coldest part of the cores. Einally, the filament model may not 
be a particularly good description of isolated clumps (such as 
LDN 183) and this may lead overestimation of the line-of-sight 
temperature gradients. Eurther analysis of p variations would re¬ 
quire detailed modelling of all the individual clumps. 

5. Conclusions 

We have used Herschel and Planck observations of dust emis¬ 
sion to examine the dust opacity spectral index in a number of 



Fig. 27. Temperature and spectral index values towards the 
submi llimetre clumps catalogued in |Montillaud et al.| ( |2015, in| 
press] ). The blue and red contours correspond to starless and 
protostellar clumps respectively. There are five equally spaced 
contours between 15% and 90% of the peak density of points. 
The background grayscale image and the colour bar show the 
distribution of unclassified clumps (logarithm of the number of 
clumps per bins of AT = 0.14K and isp = 0.014). The white 
cross shows the peak of this distribution. The plotted data have 
been smoothed with a Gaussian with EWHM equal to 1.0 K and 
0.1 units in p. 


fields containing dense clumps. The fields, typically ~ 40' in 
size, were originally selected based on the signature of cold dust 
emission but they also cover some more diffuse areas with higher 
dust temperatures. Our goal was to measure the variations of the 
spectral index and to look for correlations with environmental 
factors. We have presented maps of T and p derived with differ¬ 
ent combinations of Herschel, IRAS, and Planck data. The study 
led to the following conclusions: 

- The observed dust opacity spectral index shows clear spatial 
variations that are anti-correlated with temperature and cor¬ 
related with column density. Based on IRAS and Planck ob¬ 
servations the median values of field averages are T = 16.1 K 
andySpiR = 1-84. Using Herschel data, we find in many indi¬ 
vidual cold clumps values up to yS ~ 2.2 (v < 353 GHz). 

- The millimetre emission exhibits a clear excess relative to 
the modified blackbody fits that describe the submillimetre 
data. Two component fits to the surface brightness data aver¬ 
aged over all fields yield values ySpiR = 1.91 and p^ = 1.66. 
As expected for a sample of high column density fields, ySpiR 
is much higher than values reported for diffuse medium. On 
the other hand, y^mm is close to the values found in Planck 
studies of diffuse regions at high Galactic latitudes. 

- The amount of CO emission in the Planck bands is one of 
the main sources of uncertainty for the y^mm values. The issue 
could be resolved by directly mapping the CO(2-l) emission 
with ground-based radio telescopes. 

- The spectral index values are lower close to the Galactic 
plane, possibly because of the increased mixing of dust at 
different temperatures. We do not observe any significant de¬ 
pendence on Galactic longitude. 

- We have examined how p estimates are affected by various 
error sources and the data sets used. We conclude that our re¬ 
sults on p variations are robust concerning the detected local 
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spatial variations (vs. dust temperature and column density) 
and the variations as a function of wavelength. 

- Because temperature gradients in optically thick sources lead 
us to underestimate the intrinsic p, the true variations of the 
spectral index and hence of the dust properties may be even 
more pronounced. Quantitative estimates of intrinsic p will 
require detailed radiative transfer modelling of the observa¬ 
tions. 
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Notes. The aliases refer the followin g catalogues: Arch = [Desert et al.|{|2008), B = |Barnar dl(M9l;|Barnard et al.|fT927l, BDN = |Bemes |(T977l, CB =|Clemens & Barvainis||T988l, FeSr = 

IFeitzinger & StueweH1984) HMST = [Hartley et al.tfl^ , LDN = |LyndsH1962) , MBM = |Magnani et al.H1985], MLB^lMyers et al.MigW SLDN =|Sandqvist & LindroosH1976] . Most entries 

are looked up from |Dutra & Bica| i |2002] . Names starting with PCC and the entries Sl-SlO refer to names used in |Juvela et al. H2010] and |Planck Collahoration XXII pOl 1] forsome of the Herschel 
fields. 
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Fig.A.l. Cross sections of the cylindrical radiative transfer models. The peak optical depths are ry = 10.0 in the left and ry = 50.0 
in the right panel. Each frame shows the Ty profile, the actual dust temperature in the model midplane (Tsd), and the colour 
temperature Tq and spectral index p derived from the simulated surface brightness maps. 


Appendix A: Radiative transfer modeis 

Modelling was used to assess the possible impact that radiative transfer effects have on the observed spectral index values. The 
model consists of a cylindrical filament with a Plummer-like radial density profile 


n{r) = 


riQ 


\ + {rlRcr 


(A.1) 


where r is the distance from the symmetry axis, hq is the central density, and Rq determines the extent of the inner flatter part of 
the density profile . The filament is illuminated externally by interstellar radiation field ( [Mathis et al.|1983| ). We use the jOssenkopfj 
& Henning (19941 dust model for coagulated grains with thin ice mantles accreted over 10^ years in a density of ^ = 10^ cm“'^. To 


remove any ambiguity from the comparison of the observed and the intrinsic p values, the model is modified so that p is equal to 

1.9 at all wavelengths X > lOOyum. Note that this is the intrinsic value of the dust opacity spectral index. Because of line-of-sight 
temperature mixing, th e apparent value estimated from surface brightness data can be expected to be lower ( jShetty et al.||2009b[ 
[Juvela & Ysard|2012bj ). 

The filament is observed in a direction perpendicular to its symmetry axis. The radiative transfer calculations are used to solve 
the three-dimensional temperature structure of the filament and to calculate surface brightness maps at 160, 250, 350, and 500yum. 
The surface brightness data are used to derive maps of observed colour temperature and spectral index. In the calculation we use 
10% and 2% relative uncertainties for the PACS channel and the SPIRE channels, the same as in the case of actual observations. The 
maps consist of 256x256 pixels corresponding to the 256^ cell discretisation of the 3D model cloud. The models are characterised 
by their total optical depth in V band that, measured through the filament in the perpendicular direction, is either 10.0 or 50.0. To 
make the model more concrete, we assume for the models a linear size of 10 pc and a distance of 400 pc. With this scaling, the 
parameter Rq is equal to 0.1 pc and is determined by the selecte d opti cal depths. Each pixel corresponds to ~20" and the final 
maps are convolved to a resolution of 40.0" before analysis. Eigure [AT] shows cross sections of the modelled filament in column 
density (V band optical depth), dust temperature, and the colour temperature and the spectral index estimated from the surface 
brightness data. 

Eor the riV) = 10 model, the apparent spectral index decreases to 1.75, 0.15 units below the intrinsic p value of 1.9. Eor the 
riV) = 50 model, the minimum observed value is only 1.10. This demonstrates the difficulty of making a direct connection between 
the apparent spectral indices and the intrinsic dust properties in the case of optically thick clouds. However, because the line-of- 
sight temperature variations always imply a decrease in the apparent spectral index, this increases the significance of the observed 
correlation between column density and spectral index. At least in some regions along the line-of-sight, the actual dust opacity 
spectral index should be much higher than the observed values. 

We calculated two additional cases where the dust properties were modified in the central regions with density exceeding 25% 
of its maximum value. The modification is applied only to wavelengths T > 40 yum and it is determined by two criteria. Eirst, the 
dust opacity is set to twice its original value at 250yum and, secondly, the long wavelength spectral index p' is set to a value of either 
1.7 or 2.2. These are ad hoc models where the dust opacity has discontinuity at 40yum. This has little practical importance because 
grains absorb energy at wavelengths shorter than this limit and correspondingly emit the energy at much longer wavelengths. 

The first result is that the observed (i.e., values derived from surface brightness data) p values towards the filament centre 
are lower in both cases of modified dust. We interpret this as an indication that the increased line-of-sight temperature variations 
associated with the enhanced submillimetre opacity outweigh the effect of AyS -0.2 changes of the intrinsic spectral index. Eor the 
Tv = 10 model, the minimum p values were 1.41 and 1.67 for the p' = 1.7 and p' = 2.2 cases, respectively. Both are lower than the 
previous apparent value of 1.75. Eor the ry = 50 model, the corresponding minimum values are 0.57 and 0.67, much lower than the 

1.10 obtained before modifications to the dust model. 

Eigure |A^ shows selected SEDs from the Ty = 50 models described above. The uppermost three SEDs are observed towards 
the filament centre for the three different dust models. The two lower SEDs are obtained at different radial offsets, from a region 
that is unaffected by the P' changes that were applied to the inner part of the filament. We use arrows to indicate the difference 
between the fitted modified blackbody curves and the actual data points. Towards the optically thick centre, the data do not follow 
a single modified blackbody spectrum. The errors are largest at 160yum, because of its smaller weight in the fit (assuming 10% 
uncertainty instead of the 2% at the longer wavelengths). However, the fits also systematically overestimate the 250yum intensity 
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Fig.A.2. Selected SEDs from ry = 50.0 models. The uppermost three spectra are towards the centre of the filament, using the 
original dust model (black line) or using 13=1.1 dust (blue line) or JS = 2.2 dust (red line) in the central part of the filament. The 
lower spectra (magenta and green) correspond to larger offsets (4.4 and 8.7 arcmin) from the symmetry axis that are not directly 
affected by the dust modifications. The differences between the measured points and the modified blackbody curves are indicated 
with vertical arrows. The length of the arrow is ten times the difference between the data point and the fitted modified blackbody 
curve. 


and underestimate the 350yum intensity, this being related to the low values of the fitted spectral index. The SEDs for both modified 
dust cases are clearly colder than the original SED because of the higher opacity at wavelengths where the dust grains emit most of 
their energy. 

The spectral index was modified at T > 40//m. In the 100-500 jim range, the [3 difference corresponds to a factor of two difference 
of relative opacity. Nevertheless, the observed spectra for the {3' = 1.7 and {3' = 2.2 cases are practically on top of each other over 
the whole wavelength range probed by Herschel observations. After the opacity increase, the central part of the filament is very 
cold with dust temperatures close to 5 K. As a result, even at 500yum the emission per unit dust mass is more than a factor of ten 
lower than for the ~17 K dust on the cloud surface. This means that the {3 changes at the cloud centre are masked by the warmer, 
more strongly emitting envelope. In the model the difference of intrinsic (3 values in the filament centre becomes evident only at 
mm wavelengths. In the Herschel range, the spectral index determination is clearly very sensitive to temperature variations. If the 
160yum data point were omitted from the fits, the spectral index values of 1.16, 0.52, and 0.62 would increase to 1.37, 0.75, and 0.93 
respectively, the three values corresponding to the three models: one with original dust, one with (3' = l.l dust in the centre, and one 
with p' = 2.2 dust in the centre. The difference between the observed and intrinsic p (for given column density) would also decrease 
if the column density peak were less steep or if the observations were averaged over a larger beam. At lower column densities, at 
larger offsets from the symmetry axis of the cylindrical filament, the dust SED is well-described by a single modified blackbody 
curve, with the observed p approaching the intrinsic p value. 

In spite of such a strong effect in the above model, the observations rarely show any decrease of p towards the dense regions. 
A weak anti-correlation between p and column density was seen only in a few fields. These incl ude G82.65-2.00 w here, based on 
the observed dust emission, the column density exceeds A^(H 2 ) = 2 x 10^^ cm“^ in many places (Juvela et al. 2015). WISE 12yum 
data show clear absorption along the full filament, also suggesting that the visual extinction is at least Ay ~ 20 mag. The average 
properties of the filament would therefore be between the two models above but on average closer to the ry = 10.0 case. The large 
scale decrease of p towards the filament is at the level of AyS = 0.1 and could easily be explained by temperature variations within the 
beam. However, even in G82.65-2.00 the T -yS relation shows an anti-correlation at the smallest scales towards the densest clumps, 
contrary to the behaviour in the models. 

The temperature mixing in regions like G82.65-2.00 is probably very complex compared to the simple models above and the 
net effect on the observed spectral index can be qualitatively different. The models above showed a maximal effect because we 
compared emission between the cold centre and the fully illuminated cloud surface. In real clouds the dense regions are shielded by 
diffuse envelopes whose extinction decreases temperature contrasts at small scales. The drop in the observed p would be smaller if 
the column density profile were less steep or if the data were smoothed by a larger beam. If the emission region is inhomogeneous, 
elongated, or consists of completely separate clumps along the line-of-sight, the radiative tran sfer effects would be reduced. Internal 
heating of some clumps is not able to explain the general anti-correlation between T and p ( [Malinen et al.|20lT| [Juvela & Ysard 


|2012b| ). Therefore, it seems that we are observing real changes in the intrinsic dust properties that, remarkably, are strong enough 

to be seen above the opposite effects caused by temperature variations. _ 

The models can also be used to examine the dependence on the wavelengths used. Eigure [A3] shows temperature and spectral 
index profiles of the riy) = 10 model that are estimated using five different wavelength combinations. As expected, p values are 
least biased when long wavelengths, 250-850yum, are used. However, in the case of this model, even the combination of PACS and 
SPIRE bands 160-500yum results in values that are lower by no more than AyS = 0.07. This suggests that over most of the observed 
fields the bias resulting from temperature variations is not significantly different for the different combinations of frequency bands. 
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Fig. A.3. Temperature and spectral index profiles of the t(V) = 10 model observed using different wavelength combinations. In 
order of decreasing colour temperature these are: 100, 350, 550, and 850yum (black curves), 160, 250, 350, and 500yum (green 
curves), 160, 350, 550, and 850yum (blue curves), 250, 350, and 500yum (cyan curves), 250, 350, 550, and 850yum (red curves). 
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Fig.B.l. Modified blackbody fits to Planck and IRIS data in fields G0.02+18.02 and G0.49+11.38. The frames a, c, and d show 
250yum surface brightness, colour temperature, and spectral index fitted to 3000 GHz (JRAS 100yum), 857 GHz, 545 GHz, and 
353 GHz data, respectively. Frame b shows the higher resolution Herschel t( 250yum) map. Frame e is the surface brightness residual 
between the 217 GHz observation and the above fit. The last frame shows the difference between spectral indices derived with fits to 
3000 GHz-353 GHz and to 3000 GHz-217 GHz data. The fields in this and the following figures are in order of increasing galactic 
longitude. 


Appendix B: Fits of T, p with IRAS and Planck data 

The figures [bT| - [B. 59 1 show results of modified blackbody fits of the 100yum IRAS data and the Planck 857 GHz-353 GHz observa¬ 
tions. The bottom frames show the 217 GHz residuals (observation minus the prediction of the fitted model) and the difference in 
spectral index when the lowest included frequency is 353 GHz or 217 GHz. For the field G300.86-9.00 the plot was already shown 

in Fig. 12 
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Fig.B.2. Continued... Fields GO.02+18.02 and GO.49+11.38 
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Fig.B.3. Continued... Fields Gl.94+6.07 and G2.83+21.91 
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Fig.B.4. Continued... Fields G3.08+9.38 and G3.72+21.02 


I 

s 


I 

i 



Fig.B.5. Continued... Fields G4.18+35.79 and G6.03+36.73 
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Fig.B.6. Continued... Fields G9.45+18.85 and GlO.20+2.39 



Fig.B.7. Continued... Fields G20.72+7.07 and G21.26+12.11 
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Fig.B.8. Continued... Fields G24.40+4.68 and G25.86+6.22 




Fig.B.9. Continued... Fields G26.34+8.65 and G37.49+3.03 
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Fig.B.lO. Continued... Fields G37.91+2.18 and G39.65+1.75 
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Fig.B.ll. Continued... Fields G62.16-2.92 and G69.57-1.74 



Fig.B.12. Continued... Fields G70.10-1.69 and G71.27-11.32 
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Fig.B.13. Continued... Fields G82.65-2.00 and G86.97-4.06 
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Fig.B.14. Continued... Fields G89.65-7.02 and G91.09-39.46 
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Fig.B.15. Continued... Fields G92.04+3.93 and G92.63-10.43 
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Fig.B.16. Continued... Fields G93.21+9.55 and G94.15+6.50 
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Fig.B.lS. Continued... Fields G105.57+10.39 and G107.20+5.52 
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Fig.B.19. Continued... Fields G108.28+16.68 and G109.18-37.59 






o 

S 



o 

I 

e 





208 
200 
19 2 
184 

•X 

176 
16 8 
160 

15 2 
0100 
0075 
0050 \ 
0 025 £ 

0000 S 

-0 025 I 
-0 050 
-0075 


Fig.B.20. Continued... Fields G109.80 -f 2.70 and G110.62-12.49 
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Fig.B.21. Continued... Fields GllO.89-2.78 and GllO.80+14.16 
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Fig.B.22. Continued... Fields Gil 1.41-2.95 and G115.93-^9.47 
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Fig.B.23. Continued... Fields G116.08-2.40 and G126.24-5.52 
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Fig.B.24. Continued... Fields G126.63 -f 24.55 and G127.79 -f2.66 
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Fig.B.25. Continued... Fields G128.78-69.46 and G130.37+11.26 



Fig.B.26. Continued... Fields G130.42-47.07 and G131.65+9.75 


42 





















































M. Juvela et al.: Galactic cold cores VI. Dust opacity spectral index 




Fig.B.27. Continued... Fields G132.12+8.95 and G139.60-3.06 
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Fig.B.28. Continued... Fields G141.25+34.37 and G149.67+3.56 


43 





























































M. Juvela et al.: Galactic cold cores VI. Dust opacity spectral index 



Fig.B.29. Continued... Fields G150.47+3.93 and G151.45+3.95 
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Fig.B.30. Continued... Fields G154.08+5.23 and G155.80-14.24 
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Fig.B.31. Continued... Fields G157.08-8.68 and G157.92-2.28 
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Fig.B.32. Continued... Fields G159.23-34.51 and G159.12-14.30 
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Fig.B.33. Continued... Fields G159.34+11.21 and G161.55-9.30 
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Fig.B.34. Continued... Fields G163.82-8.44 and G164.71-5.64 
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Fig.B.35. Continued... Fields G167.20-8.69 and G168.85-10.19 



Fig.B.36. Continued... Fields G171.35-38.28 and G173.43-5.44 
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Fig.B.38. Continued... Fields G181.84-18.46 and G188.24-12.97 
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Fig.B.39. Continued... Fields G189.51-10.41 and G195.74-2.29 



Fig.B.40. Continued... Fields G198.58-9.10 and G202.23-3.38 
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Fig.B.41. Continued... Fields G202.02+2.85 and G203.42-8.29 
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Fig.B.42. Continued... Fields G205.06-6.04 and G206.33-25.94 
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Fig.B.43. Continued... Fields G210.90-36.55 and G212.07-15.21 



Fig.B.44. Continued... Fields G215.37-3.04 and G215.44-16.38 
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Fig.B.45. Continued... Fields G216.76-2.58 and G218.06+2.12 
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Fig.B.46. Continued... Fields G219.36-9.71 and G219.29-9.25 
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Fig.B.47. Continued... Fields G227.95-2.98 and G247.55-12.27 



Fig.B.48. Continued... Fields G253.71+1.93 and G255.33-4.88 
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Fig.B.49. Continued... Fields G258.90-4.10 and G265.04+6.08 
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Fig.B.50. Continued... Fields G265.60-5.82 and G268.21+2.02 
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Fig.B.51. Continued... Fields G271.06+4.84 and G271.51+5.14 
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Fig.B.52. Continued... Fields G276.78-i-l.75 and G298.31-13.05 
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Fig.B.53. Continued... Fields G299.57+5.61 and G300.61-3.13 



Fig.B.54. Continued... Fields G300.86-9.00 and G315.88-21.44 
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Fig.B.55. Continued... Fields G320.84+5.09 and G325.54+5.82 
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Fig.B.56. Continued... Fields G332.70+6.77 and G334.65+2.67 
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Fig.B.57. Continued... Fields G339.22-6.02 and G341.18+6.51 
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Fig.B.58. Continued... Fields G343.64-2.31 and G344.77+7.58 
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Fig.B.59. Continued... Fields G345.39-3.97 and G358.96+36.75 
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Fig. C.l. Temperature and spectral index fits of field G0.02+18.02 using Herschel data. Frames are: colour temperature and spectral 
index with SPIRE data and;^^ minimisation (frames a - h), 250)um optical depth and colour temperature for SPIRE data and yS =1.8 
(frames c-d), colour temperature and spectral index with SPIRE data and MCMC calculations (frame e-f) and the corresponding 
error maps (frames g-h), MCMC results for 160-500yum fits (frames i-j) with corresponding error estimates (frames k-l). The last 
frames show the effect of zero point uncertainty in the three-band fits (frames m-n) and four-band fits (frame o-p) ix^ fits) based on 
Monte Carlo simulation usingfits. Temperature maps have contours drawn at intervals of 2.0 K, starting at 10.0 K. Spectral index 
maps have contours at intervals of 0.2, starting at 1.0. 


Appendix C: Maps of (r,yS) with Herschel data 

The figures |C.lf|C.ll?] show the colour temperature and spectral index maps calculated with Herschel data. The figures include 
results derived with both 160-500 yum and 250-500yum data. 
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Fig. C.2. Continued... Field G0.49+11.38 
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Fig. C.3. Continued... Field Gl.94+6.07 
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Fig. C.4. Continued... Field G2.83+21.91 
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Fig. C.5. Continued... Field G3.08+9.38 
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Fig. C.6. Continued... Field G3.72+21.02 
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Fig. C.7. Continued... Field G4.18+35.79 


63 


TtJ=a.S||K> ^ «4^, hUj TtJ=I.S|iK^ 























































































M. Juvela et al.: Galactic cold cores VI. Dust opacity spectral index 



o 

I 


JiSb- 

PAtj2iHS0i 



aSEB*- JLiSb-' JJll-4' 



l-Qp 

?.i 

IS 

43 

3.S 

la: 

21 

14 
4gi 

25 

25 

2D 

IQ 

15 
li 
3 4 
1J 
10 
OJ. 


2SB e<- J3a €* 23a 4* 


6 

TZ 


€ 




k 



23BE+ 33a 4" 

PAg2(KMi 


Fig. C.8. Continued... Field G6.03+36.73 
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Fig. C.9. Continued... Field G9.45+18.85 
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Fig. C.IO. Continued... Field GlO.20+2.39 
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Fig. C.15. Continued... Field G26.34+8.65 
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Fig. C.18. Continued... Field G39.65+1.75 


Fig. C.19. Continued... Field G62.16-2.92 
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Fig. C.20. Continued... Field G69.57-1.74 
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Fig. C.21. Continued... Field G70.10-1.69 
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Fig. C.23. Continued... Field G82.65-2.00 
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Fig. C.32. Continued... Field G98.00+8.75 
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Fig. C.33. Continued... Field G105.57+10.39 
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J"' 'Uj TtJ=!.S||K) JJ" in, hUj TCJ = i.SHK:- 
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TXJ=j.S||K> di, hUj T(.J = i.SMK.> 
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Fig. C.36. Continued... Field G109.18-37.59 
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TtJ = i.Sl|K> Jr' >1^. tcj=i.smk:- 
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Fig. C.39. Continued... Field G110.89-2.78 
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Tti = I.S||K> liH, hUj TCJ=1.S|iK:- 
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Ji“' '4^, ntJj T(i = a.S||K> -Jr’ di, hUj TCJ=i.SMK:i 
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Fig. C.42. Continued... Field G115.93+9.47 
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Fig. C.43. Continued... Field G116.08-2.40 
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TtJ = a.S||KJ J?" TXJ=i.SHK:- 
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Fig. C.44. Continued... Field G126.24-5.52 
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Fig. C.45. Continued... Field G126.63+24.55 
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Fig. C.46. Continued... Field G127.79+2.66 
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Fig. C.47. Continued... Field G128.78-69.46 
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TtJ = a.S||K> df" 'tJj TCJ=i.SMK:- 
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Fig. C.48. Continued... Field G130.37+11.26 
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Fig. C.49. Continued... Field G130.42-47.07 
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Jr’ hLTi TtJ = i.S||K:i JT' +Lrj ■lXJ=i.siiK:- 
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Fig. C.51. Continued... Field G132.12+8.95 
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TCJ=a.S||K> +Lrj TCJ=i.SUK:- 
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Fig. C.52. Continued... Field G139.60-3.06 
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j'" ^LTi TtJ = a.S||K:i ilH, hUj TCJ = i.SMK:i 
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Fig. C.54. Continued... Field G149.67+3.56 



16 5 

I«4 

H& 

U4] 

Mi 

Mfl 




&36 


0.22 


OS I 


185 

18D 

in 

1 741 
165 


H B+ 64. B'- 

itAg?(K»k 


til 2» 

I +S6 1- 

4-54.B* 

+56.3* 

_ s* 

I -1^11* 

+54 B* 
+^13* 

1-56.2* 

I 1’ 

1 

^ +54 B* 
f54.fl.* 

-1-56.3* 
+K 2* 
I 4.S^1* 

* 

+54 B* 



66.4* 66fl* 

PAgaoook 





66 4* 6641' 

RAg?CH»h 



66 4* 66 6' 

M.g2(K»^ 


2.2 

2.0 

16 

14 

IJ 

1 £Q 

162 

1T5 

I6t - 

161 ^ 

164 

14f 



Fig. C.55. Continued... Field G150.47+3.93 
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TXJ = i.S||K:' ^ ilH, Mi TCJ^i.SliKJ 


















































































































M. Juvela et al.: Galactic cold cores VI. Dust opacity spectral index 


I ^ 3* 





+S4.1» 


I t^S4.3* 

Li 

Y t-M.a* 

I 

*54 1+ 

I 

^1:6? 
g +5J.1* 


*54.D* 



6? 6' €T i* 

PAtjiaCfQi 



t? t" €7 a* 
RAga(H»h 


L. 



^3 




m 

1 

1 » 


T^KQijmil 




.1“' 


t? 6' 67 a» 


35 

3.D 

ss 

2a 

15 

IQ 

44 

ta 

IS 

J2 

25 

14 

1#D 

1611 


1S6 

J44 

iia 

iJfl 

i^q 

190 
1 7& 
1 76 
16E! 
166 



6? 6' 67 a» 

KAgacKXJk 


Fig. C.56. Continued... Field G151.45+3.95 
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Fig. C.57. Continued... Field G154.08+5.23 
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Fig. C.58. Continued... Field G155.80-14.24 
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Fig. C.59. Continued... Field G157.08-8.68 
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j“' dH, Tti = I.S||K:. JT' TCJ = i.S|iK:- 
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Fig. C.60. Continued... Field G157.92-2.28 
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Fig. C.61. Continued... Field G159.23-34.51 
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Fig. C.62. Continued... Field G159.12-14.30 
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Fig. C.63. Continued... Field G159.34+11.21 
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Fig. C.64. Continued... Field G161.55-9.30 
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Fig. C.65. Continued... Field G163.82-8.44 
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Fig. C.66. Continued... Field G164.71-5.64 
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Fig. C.67. Continued... Field G167.20-8.69 
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hLTj TCJ = a.S||K:i Jr“ *Lrj TCJ=i.SMK:- 
















































































































M. Juvela et al.: Galactic cold cores VI. Dust opacity spectral index 



I m r 

^ ... 

I 

^ +31 
... 
+JS 6* 

|... 

^ +JS f}* 

I ... 

+-31.B.* 

g.i. 

= 1-31.B.’ 

I 

* i-il.7* 

+^33 6* 


US 
Ul ? 
M- a ^ 
!£■ 1 
15-4 

M 

2.0 

4 

IB k 

1 7 


17 B 
i«a 

15-4 rT 

,J^ 

14 T H 
14 B 

20 

14 


e4 



*« 

RAg}(H»h 


21 

2.0 


li 

ISQ 





£4 2^ 


0.^ 

0.4 ^ 

■n 

0.? 2 

OS I 

o. ^ 

■n 

0.4 

40 

44 

4.0 

AO - 
IS ^ 
24 

24 

24 

23 

la 

J€ 

14 

Its 

1 7Ei 
Itfl , 

isi 

3.47 



€4 1 * 
PAgjiKMk 


17 B 
16 5 
16 B 

H5 

U4 

14 E 

042 

041 

0.40 

Ol? 

0^ 

024 

0 2? 

0 22 
3^5 
]« 
laE 

ISO 

175 

17B 

165 


Fig. C.68. Continued... Field G168.85-10.19 
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Fig. C.69. Continued... Field G171.35-38.28 
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Fig. C.70. Continued... Field G173.43-5.44 
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Fig. C.71. Continued... Field G174.22+2.58 
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Fig. C.72. Continued... Field G176.27-2.09 


IrlT 2* 

4-1 T.l* 

8 

w ^16.a+ 

I 

4-iie* 

*li7* 
4-1 IS' 

4-17.1* 

5 4-llr.(r 

1 

-I-I6.B.* 

4-l*.7* 

-1-17.2* 

■kITI* 

I *17 &* 

I" “ 

*li 

+ I9.7* 
4-1 T 2* 

4-17.1+ 

I 4-IIfr* 

4-li.l+ 

-l-li.7* 





14 I ,’7’ 
14 41 ^ 


12 1 
194 
17 € 
199 

199 ' 

I 

152 
14 4 

ll€ 
10 2 


E 

1*9 

15? 

217 
2.10 
2 9i - 
IK 

192 


1« ^ 

I2i b 

i«i| 

oii g 

t: 

050 

50 

35 _ 

iic 

30 “ 

?7 ^ 

24 

21 

250 

725 

2» g 

17& ^ 


11^ 

lOfi 

19Q 

192 

1 7E I 

i*e 

lil 

154 


17 0 

U4 

IS B 
159 
14^ 
0.3* 


0.2JI 


0.22 


7^5 
IK 
19S 
192 1^-. 

m 

i*fl 

191 


71.2+ 71.0+ 70 9* 

[IAg?EH»k 


71.2+ 71.0+ 70 9* 

RAg?Di»h 


71.2+ 71.0+ 70 9* 

[U.g?llOO| 


71 2+ 71 0+ 70 9* 

PAg?Dogi 


Fig. C.73. Continued... Field G181.84-18.46 


96 


TtJ=a.S||K> JT' +tJj TCJ^i.SMK) 

























































































































M. Juvela et al.: Galactic cold cores VI. Dust opacity spectral index 


I 

+H 6* 

i-T6.a» 

I 

♦ It 2* 

I 

I 

■hl-i.ri.- 



»€■ 7a-2' 




0 9- r- 

•n 

0? % 
ot I 

05 g 

TI 

Q.4 

4 5 
42 
35 

34 ^ 

53- C 

34 

i.1 

i2t 
?L4 
?4* M 

190 

ItB 

195 

lt4 

1 

170 ^7 

1«& 

IM 


»e* 79 2" 


Tit* 

RAgaiii»h 


23 

2,2 

21 .-, 
^4 V 

If 

22 

i\ 

2> 

li -r7 

19 


1 TO 



h ■ 

1 r 


^1 


■ 


m 

4Jj ^ 

_ 



17 5 
17.0 ^ 
16 5 5 

II 

160 ^ 
IS 5 
ISO 
5 39 

O lT ^ 

5.36 



Tit* 792* 

PAg?(K»k 


195 

190 

L 

190 
1 75 

1 To 


Fig. C.74. Continued... Field G188.24-12.97 
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Fig. C.75. Continued... Field G189.51-10.41 
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Fig. C.76. Continued... Field G195.74-2.29 


+oa.&+ 
*oa 5* 
g +«4* 

= *oa3+ 

i +o?2* 

+M1* 

fO4.0* 

^ fM.i- 
§ toa.^- 
3 +M 3* 
JS -HM.i* 
tOfl.l- 

*oa [>+ 

tOfl.t- 

*oa.5+ 

S 

^ 1-M.3- 

S *oa. 2 + 
+oa 1 + 
+M0^ 
*oa &+ 
_ +«5* 
§ -HM ^+ 

3 *06 3+ 

I +«s* 

firf 1+ 

toa.D* 





^4^ tfsti- 


PAgaaook 



£a4* 4f541' 

RAg?CK>0h 



a&4> 


(ugauoo^ 





SB 4* 

PAgJOOQI 


Fig. C.77. Continued... Field G198.58-9.10 
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Fig. C.78. Continued... Field G202.23-3.38 
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Fig. C.79. Continued... Field G202.02+2.85 
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Fig. C.80. Continued... Field G203.42-8.29 
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Fig. C.81. Continued... Field G205.06-6.04 
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Fig. C.82. Continued... Field G206.33-25.94 
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Fig. C.83. Continued... Field G210.90-36.55 
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Fig. C.84. Continued... Field G212.07-15.21 



Fig. C.85. Continued... Field G215.37-3.04 
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Fig. C.86. Continued... Field G215.44-16.38 
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Fig. C.87. Continued... Field G216.76-2.58 
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Fig. C.88. Continued... Field G218.06+2.12 
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Fig. C.89. Continued... Field G219.36-9.71 
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Fig. C.90. Continued... Field G219.29-9.25 
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Fig. C.91. Continued... Field G227.95-2.98 
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Fig. C.92. Continued... Field G247.55-12.27 
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Fig. C.93. Continued... Field G253.71+1.93 
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Fig. C.94. Continued... Field G255.33-4.88 



Fig. C.95. Continued... Field G258.90-4.10 
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Fig. C.96. Continued... Field G265.04+6.08 
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Fig. C.97. Continued... Field G265.60-5.82 
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Fig. C.98. Continued... Field G268.21+2.02 
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Fig. C.99. Continued... Field G271.06+4.84 
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Fig. C.IOO. Continued... Field G271.51+5.14 
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Fig. C.lOl. Continued... Field G276.78+1.75 
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Fig. C.102. Continued... Field G298.31-13.05 
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Fig. C.103. Continued... Field G299.57+5.61 
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Fig. C.104. Continued... Field G300.61-3.13 
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Fig. C.105. Continued... Field G300.86-9.00 
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Fig. C.107. Continued... Field G320.84+5.09 
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Fig. C.108. Continued... Field G325.54+5.82 
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Fig. C.109. Continued... Field G332.70+6.77 
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Fig. C.llO. Continued... Field G334.65+2.67 
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Fig. C.lll. Continued... Field G339.22-6.02 
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Fig. C.112. Continued... Field G341.18+6.51 



Fig. C.113. Continued... Field G343.64-2.31 
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Fig. C.114. Continued... Field G344.77+7.58 
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Fig. C.115. Continued... Field G345.39-3.97 
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Fig. C.116. Continued... Field G358.96+36.75 
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Fig.D.l. Modified blackbody fits in field G0.02+18.02 using the combination of Herschel and Planck data. The uppermost frames 
show the fitted intensity at 250yum, and the colour temperature and spectral index maps. The relative residuals (observation minus 
model, divided by model prediction) are shown in frames d and / for the Planck bands of 353 GHz and 217 GHz. Frame / shows 
the median values of the residuals for the three Herschel bands (red symbols) and the four Planck bands (blue symbols). The error 
bars correspond to the median error estimate of the surface brightness over the map (mainly the assumed calibration errors and the 
uncertainties of the CO corrections). The lower data points at 850yum and 1380yum {Planck bands 353 GHz and 217 GHz; cyan 
symbols) correspond to twice the default CO correction. 
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Fig.D.2. Continued... Field G0.49+11.38 


Appendix D: Combined {T,/3) fits to Planck and Herschel data 

The joint fit of Herschel (250 yu m-5Q 0 /im) and Planck (857 GHz-217 GHz) data was carried out only for the 53 fields with estimated 
distances d <<400 pc (see Sect. 1^. For comparison, the figures also include 353 GHz and 217 GHz residuals from fits where we 
assume twice the default CO correction (see Sect. |2.2| ). The larger correction could in some cases eliminate or even reverse the sign 
of the 217 GHz residual. However, for example in the case of field G6.03-r36.73, grou nd-based observations confirm that the line 
ratios assumed in the default CO correction are of the correct magnitude (see Sect. |2.2| ). 
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Fig.D.3. Continued... Field Gl.94+6.07 
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Fig.D.4. Continued... Field G2.83+21.91 
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Fig.D.5. Continued... Field G3.08+9.38 



Fig.D.6. Continued... Field G3.72+21.02 
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Fig.D.7. Continued... Field G4.18+35.79 
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Fig.D.8. Continued... Field G6.03+36.73 
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Fig.D.9. Continued... Field G9.45+18.85 
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Fig.D.lO. Continued... Field G20.72+7.07 
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Fig.D.ll. Continued... Field G21.26+12.11 
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Fig.D.12. Continued... Field G24.40+4.68 
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Fig. D.13. Continued... Field G25.86+6.22 
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Fig.D.14. Continued... Field G26.34+8.65 


125 





























































M. Juvela et al.: Galactic cold cores VI. Dust opacity spectral index 



+^72.?’ 
+72 fl' 


GIOS 2S+LS 6S 


M T 



44 


M 

7J 

S 

2* 

» 

4 

* 

4 

-2 

-M 

*4 


k 

-ki 



I 

i 


^k4.S' 314.0- 31M+ 3k 3k4.0^ 
■U UZWHI- 




]7.0 

kO.3 

k4.0 

l^.^ - 
K 

IJ.O I- 

34.0 

]J.O 

13.3 

]3.D 

«□ 


31*4^ 3Sfe0r^ 3k 3.3- 314.4- 

lU. > 



L33 

^k4> 

2»^ 

700 

lU 

LU 

L34 



Fig.D.15. Continued... Field G108.28+16.68 
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Fig.D.16. Continued... Field G110.80+14.16 
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Fig.D.17. Continued... Field G126.63+24.55 
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Fig.D.18. Continued... Field G141.25+34.37 
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Fig.D.19. Continued... Field G149.67+3.56 
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Fig.D.20. Continued... Field G150.47+3.93 
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Fig.D.21. Continued... Field G151.45+3.95 
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Fig.D.22. Continued... Field G154.08+5.23 
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Fig.D.23. Continued... Field G155.80-14.24 
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Fig.D.24. Continued... Field G157.08-8.68 
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Fig.D.25. Continued... Field G159.23-34.51 
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Fig.D.26. Continued... Field G161.55-9.30 


131 

































































M. Juvela et al.: Galactic cold cores VI. Dust opacity spectral index 


£ 


i 


+W.a= 

+ •>4.1- 
+14.±- 

+ » r 
+>?.?■ 

+ 37 .n = 
+14.1- 

+14.4- 
*14.1- 
*14.2' 
+14. L" 
*14.0^ 
*4?.9’ 
+1M- 



94.4' 

■M 


H 

n 

n 

44 

44 

H 


-ki 


I 


■r+ 

■J. 


J4 

4^ 


^4 

9.1 


4.4 




4.1 

14 


!■» 

4.4 


\ 


-11 


-1.4 



».4‘ »4^ 

riA UI^DDOjl 


3>J 
31.6 
!i1.4 
33.2 - 
35.P I- 

34.1 
3'll 
3^ i 

34.2 


n.i 

».i 
l?.l 
K.fl i 

«.! I 

J4.Q i 

3T.l| 

33.0 

37-3i 


2 47 
7.M 
Z4t 
1« 
L41 
142 
3.44 
LHl- 




Fig.D.27. Continued... Field G164.71-5.64 
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Fig.D.28. Continued... Field G167.20-8.69 
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Fig.D.29. Continued... Field G173.43-5.44 




T 

+04.5" 
+M.4> 

+44.Z> 

+&4.1" 

+«.!■ 

+M.S- 

+ M.1- 
lllil.J' 

+P^.a' 

+i>i i," 





I 


] 

I 


RU U^AODI- 



HbA- 

TUl UL^ DDQ^ 



i.22 

Lk4 

Lk4 

T.M 

3.W “ 

3.W 

L44 

194 

L3l 



Fig.D.30. Continued... Field G203.42-8.29 
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Fig. D.31. Continued... Field G205.06-6.04 
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Fig.D.32. Continued... Field G206.33-25.94 
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Fig.D.33. Continued... Field G210.90-36.55 
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Fig.D.34. Continued... Field G212.07-15.21 
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Fig. D.35. Continued... Field G247.55-12.27 
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Fig.D.36. Continued... Field G298.31-13.05 
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Fig. D.37. Continued... Field G300.61-3.13 
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Fig.D.38. Continued... Field G300.86-9.00 
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Fig.D.39. Continued... Field G315.88-21.44 
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Fig.D.40. Continued... Field G341.18-^6.51 
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Fig. D.41. Continued... Field G344.77+7.58 
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Fig.D.42. Continued... Field G345.39-3.97 
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Fig.D.43. Continued... Field G358.96+36.75 
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